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SURGING OF AXIAL COMPRESSORS 
by 


H. PEARSON, B.A., A.F.R.Ae.S., and T. BOWMER, B.Sc. 
(Rolls-Royce Limited, Derby) 


1. INTRODUCTION 


It is fairly generally known that if, at any fixed speed of a compressor, either 
axial or centrifugal, the flow is reduced by throttling the outlet, then a point is 
eventually reached at which a complete breakdown of the air flow occurs and in 
most cases an actual flow reversal through the compressor takes place. Sometimes 
this phenoménon is gradual, taking the form of a sort of burbling, but more 
generally it takes the form of a sudden “ bang” associated with a violent shake 
of the whole foundation of the bed on which the compressor is mounted. In most 
centrifugal compressors the flow reversal which takes place stops rapidly, the 
performance recovers, and a second “ surge ” occurs if the throttling is not reduced. 
The frequency of these individual surges varies greatly according to the degree of 
throttling and many other conditions; it may be only one isolated occurrence, in 


which case it would be concluded that the working point was just, but only just 
on the surge point, or it may take the form of a rapid series of thuds indicating 
that the working point is beyond the surge. In all cases of centrifugal compressors 
known to the authors reduction of the throttling will restore the compressor to its 
normal performance. 


With axial compressors, however, reduction of throttling does not immediately 
restore the normal performance in all cases. The compressor apparently operates 
on a second “ stalled ” characteristic of reduced efficiency, and it is necessary to 
open the throttle completely and close again to restore the original performance. 
There appears to be a sort of hysteresis, and should the compressor surge when 
fitted to an engine for which control of throttling cannot be exercised, it is 
sometimes necessary to slow the engine down or possibly stop it completely in 
order to recover the normal performance. Thus while a centrifugal jet engine on 
violent acceleration will “surge” but will accelerate through this surge and pick 
up its proper performance afterwards, a similar happening on an axial compressor 
engine may necessitate a complete restart. Various schemes are known for 
avoiding such a state of affairs, such as compressor bleeding, but the best method 
of use of such schemes is to prevent surge occurring. 


Paper received May 1949 
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It is clear then that “surging” of compressors represents a definite lower 
flow limit of performance, and this limit has an unfortunate habit of occurring at 
much higher flows than is either desirable or expected. Consequently some 
knowledge of how to predict the “surge” point is desirable. To do this it is 
necessary to know why surge occurs. The following account is applicable to 
axial compressors. 

The performance of an axial compressor is the sum of the performance of 
the individual stages. Although the compressor may have been designed (although 
not necessarily so) to do the same amount of work in each stage, nevertheless as 
the flow is changed from this design point, some stages will do more work and 
others less. As the operating conditions at any stage are governed by continuity, 
i.e. by the operating conditions of those stages in front of it, it is clear that the 
overall performance of the compressor “ away from design point” is a somewhat 
complicated business involving step by step calculations. It may be stated generally 
that if pressure ratios are restricted to Jess than the design value, then as the air 
flow is reduced, the front stages tend to work at increasing incidence and therefore 
do more work, while the back stages work at decreasing incidence and do less work. 

A limit to this process is set by the stalling of the front stages; so that although 
at first the pressure ratio of the compressor may rise steeply with decreasing mass 
flow, after a time, due to the stalling of the front stages it will flatten out and 
eventually decrease again. Very roughly, when there are more stages stalled than 
unstalled this point will be reached. Thus the expected characteristics of an axial 
compressor are as shown in Fig. 1. Where the pressure is decreasing with increase 
of mass flow the slope may be said to be negative, and, working against any 
reasonable load the compressor would be stable. Where the pressure is rising 
with increasing mass flow, i.e. the slope is positive, there will evidently be a chance 
of instability depending on the characteristic of the load. In making compressor 
calculations so far it has been the custom to assume that the compressor will surge 
on the peak of its characteristic. Experimental results show that this assumption 
is well justified, and in some cases, although the form of the positive slope of the 
characteristic cannot be determined, it would appear that surging occurs at even 
higher mass flows than the peak. 


It is the object of this paper to show first that the assumption that surging will 
occur on the peak of the characteristic is not justified on simple grounds, a point 
that is probably well appreciated by those who make it; and secondly to describe 
a treatment of the problem of surging which is capable of predicting the probable 
surge point and to show that this point may well coincide with the peak of the 
characteristic or may occur at even higher flows. 


Notation 
A mean cross-sectional area of stage, sq. in. 
b constant 
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electrical capacity, farads 
potential difference, volts 
gas constant 


acceleration due to gravity, in./sec.? 


electrical current, amps. 
electrical inductance, henries 
rotational speed, r.p.m. 
pressure difference, Ib./in.? 
Pr compressor inlet pressure, Ib. /in.? 
p compressor delivery total pressure, Ib. /in.? 
Q air mass, lb. 
q electrical charge, coulombs 
R electrical resistance, ohms 
s length of stage, inches 
ye compressor inlet temperature, “K 
Ve compressor delivery total temperature, °K 
T mean stage temperature, “°K 
t time, sec. 
v air velocity, in./sec. 
V volume of the stage, cu. in. 
WwW air flow, Ib./sec. 
Z electrical impedance, ohms 
p air density, Ib./in.* 
angular velocity, rad. /sec. 
Y ratio of specific heats 


2. ELEMENTARY EXAMINATION OF STABILITY 


Consider a compressor with a characteristic at a given speed as shown in 
Fig. 1. This shows the pressure ratio p,/p, plotted against non-dimensional mass 
flow (W/T,)/p, where suffixes 1 and 2 refer to conditions at inlet and outlet 
respectively. Then suppose that this compressor is delivering air into an orifice 
whose area may be adjusted so that the working point may be at any chosen point 
of the compressor characteristic. The conditions that govern the flow through this 
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Pe 


PRESSURE RATIO 


| CHARACTERISTICS AT 
VARIOUS R.P.M/VT, 


AIR FLOW FUNCTION wWYT, 


Fig. 1 
Expected characteristics of axial compressor 


orifice are the pressure p, and temperature T, at exit from the compressor and the 
pressure at outlet from the orifice, which will be assumed to be equal to the 
pressure at inlet to the compressor p,. So to derive the correct matching of orifice 
load and compressor characteristic it is convenient to plot the compressor 
characteristic on outlet conditions. This is a familiar manceuvre in matching 
compressor characteristics and such a characteristic is shown in Fig. 2. It has 
the effect apparently of prolonging indefinitely the negative slope leg and steepening 
the positive slope leg. It is now possible to match directly to the orifice load, for 
the air flow passed may be plotted on the same diagram as a series of dotted lines 
each representing a different orifice size. The intersection of these points with 
the compressor characteristic gives the working point. 


It will be seen from Fig. 2 that, in general, there are two working points, one 
to the right of the point A and one to the left of it. The first may readily be seen 
to be stable, the second to be unstable. In the limit, testing can only be 
carried down to the point A where the two points become one, so to 
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| 


CHARACTERISTICS OF 
VARIOUS SIZED ORIFICES | 


P, 


y, / CHARACTERISTIC 
AT FIXED 


PRESSURE RATIO 


AiR Flow Function 
Fig. 2 
Matching of compressor to orifice load 


speak. If the orifice size were reduced below this, then presumably surging might 
be anticipated. 


On the other hand it may be shown that this is not likely to give a surge point in 
line with the assumption already made. If for simplicity the compressor is assumed 
to give a compression ratio greater than 1.86, then the orifice will be choked, and 
(W/T.)/p, will be constant, that is the load lines on Fig. 2 will be vertical. Thus 
instability will not now occur until the compressor characteristic itself is vertical, 
and this would then become the criterion of surge, namely that point at which 
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(W./T.)/p, is constant. As compressor characteristics plotted on inlet conditions 
are more familiar, it is preferable to translate this condition. 


If 


be constant, 


2 


then N(R) 
Pr T, 


= constant x 
T, 


This means that the surge point will be on the positive slope side of the 
characteristic whose siope is given by the above expression, i.e. the surge point is 
definitely at lower flows than the peak. In practice such calculations yield a surge 
point at flows considerably less than the peak and than the experimental results. 


3. CONSIDERATIONS OF DYNAMIC STABILITY 


It is clear that surge points at or near the peak can only be explained on the 
previous simple theory if the load line of the orifice in Fig. 2 be made much flatter. 
Under the conditions assumed, this cannot be done for steady flows, but if it is 
considered that in a surge the flow is rapidly changing, then it is important to 
consider the capacity of the system and the inertia of the air column too. If the 
compressor feeds the orifice by means of a pipe, then this will have capacity and 
inertance (or inductance to use the electrical analogy) and the effect of these will be 
to make the load line much flatter at certain frequencies. 


This explanation is fairly well known among those familiar with these 
problems, but it is unsatisfactory for two reasons. 


(i) Whatever values of capacity and inertia are used it always gives the 
surge on the positive slope of the characteristic, although maybe only just. 


(ii) It would indicate that the surge point is critical to the outlet conditions, 
i.e. whether a large, long pipe is used or not. Extensive tests seem to 
indicate that the surge point is largely unaffected by outlet arrangements. 


It is essential that if consideration is to be given to the dynamic as opposed to 
the static performance of the load, the dynamic performance of the compressor 
must be considered. In this case the air capacity and inertia of the compressor 
are by no means negligible and the compressor may be regarded as a complex 
arrangement of capacities, inertances, resistances both positive and negative, and 
so on, and it should be possible to lay down conditions for stability. The rest 
of this paper is concerned with the analysis of this problem and presents some 
representative results on a high pressure ratio compressor. 
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4. CONSIDERATION OF INTERNAL CAPACITY AND INERTANCE 


If one stage of an axial compressor is considered, we are interested in the 
relation between the different values of the pressure rise across the stage at different 
values of the mass flow, in the unsteady state, that is when the air flow is changing 
in an arbitrary manner. The extent to which the fact that the air flow is changing 
affects the performance is, to a first approximation, due to the air capacity of the 
stage and the inertia of the column of air in the stage. The problem is then 
exactly analogous to an electric circuit having capacity, inductance, and resistance. 
It appears best first to establish this analogy and afterwards to use the notation 
more familiar in electrical theory. 


5. CORRESPONDING ELECTRICAL AND AIR FLOW QUANTITIES 


Electrical Air Flow 
Potential difference (E) (volts) Pressure difference p (Ib. /in.’) 
Current 7 (amps.) Air flow W (Ib. /sec.) 
Charge q (coulombs) Air mass Q (Ib.) 
Resistance R= aE (ohms) “ Resistance ” op 
dW 
. dq “ ” dQ 
Capacity C= (farads) Capacity 
For a volume V, Q= GT 


where T=mean stage temperature 


G=Gas Constant. 


Therefore (1 ~ pat 


dp GT\ Tap 


=> GT where y=Ratio of Specific Heats. 


dp 


(henries) Inertance aw di 


dV 
Inductance L = al /dt 
If the stage is of length s and cross section A, then for a pressure difference 
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of 


Fig. 3 
Equivalent electrical circuit of a compressor stage 


dp across the length of s we have for the acceleration 


Ad 
and Q=pAs, W=pAv. 
dv 
Therefore 
_pAs 1 dW 
dt 


An elementary stage may be regarded as a circuit of the simple type shown 
in Fig. 3. If the input impedance to this circuit is regarded as negligible for the 
moment the analysis for such a simple circuit leads to a differential equation 


dl 
Lap +Ra 
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There are three possible solutions to this equation according to the actual values 
of R, L and C. 


a greater than a. In this case the circuit is critically damped and 


oscillations are not possible. 


(a) 


4L 
of frequency 


; less than -. In this case sinusoidal oscillations are obtained 


(b) le 


(ze - 


damped out gradually by the factor e-'*/?")' 
(c) =+5 =>. A Singular case, the solution being given by 


(K,t+K,) 
K, and K, being constants. 


The damping term in all these cases is the same, and so long as the resistance 
R is positive, natural oscillation, if possible, will be damped out. If R should 
become negative as it does in the case when the stage becomes stalled, then the 
“damping ” term becomes an exponential increment and there will occur rapidly 
growing oscillations or “ critical divergence ” according to whether conditions (bd) 
or (a) apply. 


The conditions in which we are interested are those that occur within the 
compressor when some stages have such a negative resistance. It remains to 
investigate the problem when there are multiple stages. 


6. EQUIVALENT CIRCUIT OF MULTI-STAGE COMPRESSOR 


For the multi-stage compressor, there is an equivalent circuit as shown on 
Fig. 4, where only the last three stages of a twelve-stage compressor are shown. 
R represents the load resistance which may be complex, i.e. consisting of resistance, 
Capacity, and inductance. C,, and L,, are the capacity and inductance of the 
twelfth stage and R,, its resistance, i.e. the slope of the pressure flow curve at the 
point of operation. The circles represent diagrammatically a generator and indicate 
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Fig. 4 


Electrical circuit analogous to a twelve stage axial compressor 


that any fluctuations of pressure from the preceding stage are multiplied by the 
compression ratio of the twelfth stage at the point of operation. All inductances 
are fixed and are directly calculable from compressor dimensions. The capacities 
are mainly fixed by the dimensions but are modified to a certain extent by the 
temperature at any stage, as outlined in Section 5. The stage resistances and 
generators all depend upon the point of operation on the stage characteristic. 


Now if a small disturbance is imagined to take place at any point in the 
compressor, this disturbance will travel either forwards or backwards through the 
stages and will either be amplified or damped out, according to the values of the 
various circuit components. The criterion of stability will be that the impedance 
at any frequency shall be greater than zero. As a criterion of surge, it may be said 
that this will occur when the impedance is zero. The minimum impedance will 
always occur at the resonant frequency and the compressor will surge when the 
resistive component at the resonant frequency is zero. 


To investigate the overall case of instability as viewed from any point within 
the compressor is a mammoth undertaking and this paper is restricted to calculating 
the impedance step by step through the compressor as viewed from the first stage. 
The impedance of the twelfth stage is calculated, viewed from the entry to the 
twelfth stage, then that of the eleventh and twelfth looking from the eleventh entry, 
and so on. When at any point the impedance becomes zero, at any frequency, 
the compressor is assumed to surge. It is appreciated that other forms of 
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instability may have been overlooked, in which, for example, the impedance 
looking towards the entry from some point in the middle of the compressor 
may be zero. 


As described in Section 5 there are two cases to consider, “ divergent 
oscillations ” and “ critical divergence.” 


(a) Divergent oscillations 


In this case since we are interested in sinusoidal oscillations use may 
be made of the complex impedance method and the impedance of a series 
circuit of resistance, capacity and inductance may be written down 
directly as: — 


| 
Z=R+ +iwL, 
where w is 2 x frequency, and i= /(-1). The impedance of the whole 
compressor may then be calculated by step by step methods, allowing 
for the fact that when passing across a stage the impedance must be 
divided by the compression ratio of the stage. 


(b) Critical divergence 


In this case the value of a small perturbation in flow J after any 
time will be 


After a short time the second term will be of no importance and we may 
write 1=K,e"". The impedance of a simple circuit containing resistance, 
capacity and inductance may then be written 


1 
Z=R+ +bL, 


i.e. the algebraical sum and not the vector sum as in the case of oscillations. 


7. METHOD OF ANALYSIS 


The overall characteristics of a twelve-stage compressor were calculated from 
assumed stage characteristics. These stage characteristics were based on a wealth 
of test data and so represented a fairly realistic compressor. The complete 
compressor characteristics so calculated are shown in Fig. 5 and the instability 
line as given by the steady flow theory, i.e. the points at which (W/T,)/p, are 
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| SURGE LINES’ 


A) CONTINUOUS FLOW THEORY 
B) OSCILLATORY FLOW 
}—C)- DIVERGING FLOW_— 


| 
| 
| 


| 


COMPRESSION RATIO R,, 


| various N | 
/ 
| 
C, | 
| 
| | | 
40 60 70 60 90 100 tO 120 130 140 


AIR FLOW FUNCTION ras 
i 
Fig. 5 


Twelve stage axial compressor characteristics 


constant is shown by the dotted line. As may be seen this criterion yields surge 
flows well to the left of the peak pressure points, which is contrary to experience. 


The impedance calculations were then made at a fixed compressor speed 
and a fixed value of the air flow for both oscillatory and divergent flow conditions, 
and a fixed value of » or b. Then the value of » or b was varied until a locus 
of the variation of impedance with frequency was obtained. If at all values this 
impedance was always positive, then the calculations were again made for a mass 
flow somewhat smaller, until the exact surge point was obtained. A typical 
result for such a calculation is shown in Fig. 6 in which resistance is plotted against 
reactance, at a fixed compressor speed and three different air flows. It will be 
seen that at (W /T,)/p,=88.5 the impedance is always positive, at (W/T,)/p,=76 
it is always negative, while at (W / T,)/p, =85.6 it changes from positive to negative. 
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Twelve stage axial compressor overall impedance at N/ / T=354 plotted on an Argand diagram 
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Nevertheless when the reactance is zero the resistance is negative and this is 
definitely a point beyond the surge. The exact surge point was found to be at 87 
in this case. 


The results for the same case assuming “ divergence,” are shown in Fig. 7. 
It will be seen that a negative resistance is obtained at (W/T,)/p,=76 but not at 
85.6 or 88.5. The surge speed on this assumption appears to be about 83.5. 


8. RESULTS OF CALCULATIONS 


The calculations for one value of mass flow and compressor speed for a 
twelve-stage compressor are very laborious, so that the magnitude of the 
investigation may be imagined. Fig. 5 summarises the results of all this work, 
and shows the loci of the surge points determined in the three ways:— 


(a) Continuous flow theory, i.e. at the points where (W/T,)/p, are constant. 
As already mentioned these points are at flows much too low to agree 
with experience. 


(5) Oscillatory flow at surge. 
(c) Diverging flow at surge. 


The most interesting point is that both conditions (b) and (c) give surge points 
very near the peak pressure points of the characteristics at high speeds and therefore 
give some justification for the assumption that the surge line coincides with the 
peak pressure locus. At lower speeds the theory of (b) and (c) seems to indicate 
that instability should occur at higher flows than the peak point, and at the lowest 
speed investigated it was not possible to find a stable point on the range in which 
calculations had been made. Therefore, if the theory on which this paper is based 
is correct, the problem would appear to be why the compressors are stable at low 
speeds at all, rather than why they stall at the peak point of the characteristic 
at high speeds. On the other hand it was found that the impedance at low speeds 
was not at all sensitive to mass flow and a small change in some of the constants 
could easily have made the compressor completely stable at all flows. 


9. DISCUSSION OF RESULTS 


The analysis of compressor stability given above, although laborious in 
practice, is still elementary in theory and before it could be regarded as a safe 
criterion of surge the following additional facts would have to be considered : — 


(a) Although account is taken of the changing flow with time, it is assumed 
that at any time the stage performance is still the same as it would be 
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Fig. 7 


3000 4000 


Twelve stage axial compressor overall impedance at N/ / T=354 


under steady flow conditions at the same value of flow. This will not 
be true, for appreciable time is taken to build up the circulation around 
the blades. The assumption is expected to be particularly wide of the 


mark for conditions under which a stage is running stalled. 


209 


| | 
| | 
| | 
| 
| | 
\ 
76 
F) | | 
0 } 
4 | | | | 
-0-1 
0 1000 2000 Pe 
| 
|_| 


H. PEARSON AND T. BOWMER 


(b) It is assumed that each stage behaves in the same way, independently of 
the manner in which previous stages are behaving. Again this is probably 
true enough when the stages in front are unstalled, but likely to be far 
from true when they are stalled. 


Both the above considerations would tend to push the surge point at low 
speeds to lower values of air flow. 


It appears to the authors that the chief value of the work described in this 
paper is to give some logical explanation of the phenomenon of surge in axial 
compressors, and to show that it is capable of predicting surge points very near 
the peak points of the characteristics, so justifying an assumption that has been 
used for a long time. Due to the complex nature of the whole process, however, 
it is doubtful if surge points can ever be predicted with any certainty. 


With regard to the “ oscillatory ” or “ critically divergent ” nature of the surge, 
although both types have been investigated, it appears to the authors that oscillatory 
instability is not very likely, due to the non-linearity of the characteristic resistance, 
which would in effect place a heavy damping on the oscillations. Oscillograph 
records taken during surging indicate no evidence of such high frequency 
oscillations, but they do show a series of sudden changes of pressure which 
gradually return to the normal value. The exact reason for this is not yet 


appreciated, and to do this would require a knowledge of the compressor 
performance while actually surging. 
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STRAIN MEASUREMENT 
BY X-RAY DIFFRACTION METHODS 


by 


G. B. GREENOUGH, M.A., Ph.D. 
(Metallurgy Department, Royal Aircraft Establishment) 


SUMMARY 


Many papers have been written on the measurement of strain by X-ray 
diffraction methods and on the interpretation of these strains in terms of stresses. 
Whereas, during the past few years, the experimental methods of determining the 
strains have remained largely unchanged, research has shown that the older 
techniques for calculating stresses from strains are not always valid. 


In this paper an attempt is made to describe some of the principles of strain 
measurement by X-ray diffraction methods to those who are unfamiliar with the 
methods. The types of stress and strain systems which may exist in polycrystalline 
metals are then considered, particular attention being paid to the effect of the elastic 
and plastic anisotropy of the individual crystals. Some indication is given as to 
how the earlier methods of interpreting X-ray strain measurements should be 
modified, but no rigid routine method is proposed for use in a general case. 


1. INTRODUCTION 


The structural engineer working in the field of aeronautics has one major 
problem, to develop a design such that each component of the structure is carrying 
a maximum load only just below that which would cause failure. At the moment 
the best that can be achieved is to calculate the maximum stress that is likely to 
occur on the component, to consult tables of data of the strengths of alloys and 
finally, to design the component such that its calculated failing load is three to 
five times the maximum to be expected. This “ safety factor” is a measure of the 
present lack of knowledge of aerodynamical stressing, of structural theory and of 
the behaviour of metals under various circumstances. The metallurgist’s contribu- 
tion to the achievement of a perfect structure is the provision of more information 
regarding the maximum stresses, either static or fluctuating, that can be applied 
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to a particular alloy and, after the component has been designed, to ensure that 
this strength is reached and is maintained throughout its life. A second contribution, 
which receives more publicity, is the development of alloys and methods of treatment 
giving greater strength than those at present in use. 


It has long been recognised that one of the factors affecting the strength of 
fabricated parts is the presence of internal stresses, that is, stresses which exist in 
the component when no external stresses are applied. Modern methods of metal 
treatment involving extrusion processes or the quenching of alloys often introduce 
these stresses adventitiously and they may have an advantageous or a deleterious 
effect on the strength of the component. On the other hand, treatments have been 
evolved for the express purpose of introducing stresses into the part to increase 
its strength; for example, the shot-peening of springs often increases their fatigue 
strength. If the effect of these treatments is to be judged quantitatively, then some 
effort must be made to measure the size and distribution of the stresses. At the 
moment the methods used depend on the measurement by one of several methods 
of the strain in the component, or of changes in strain as parts of the component 
are removed; these strains are then interpreted in terms of stresses. The whole 
field of the measurement, origin, control and effect of internal stresses has recently 
been discussed“), and one of the techniques of strain measurement there described 
was the X-ray diffraction method. 


The X-ray diffraction technique of measuring strains in metals is probably one 
of the most powerful at present employed. Many papers have been written 
describing apparatus suitabie for the purpose, but the interpretation of the measured 
strains in terms of stress systems is in an unsatisfactory state, and hitherto many 
assumptions have been made which are not always justified. In order to appreciate 
fully the significance of the results of X-ray strain measurements, it is essential to 
understand the principles involved, and it is the purpose of this paper to describe 
some of these principles to those who may be unfamiliar with X-ray diffraction 
methods. No attempt is made to discuss the details of the experimental technique 
nor the details of any mathematics involved; these have been discussed at 
length elsewhere. 


2. METAL CRYSTALS AND X-RAY DIFFRACTION 


The metals used in structural components are aggregates of crystals. Such 
crystals, or more particularly the boundaries of the crystals, can easily be seen in a 
polished and etched surface prepared on the component and examined if necessary 
under the microscope. The size of these crystals varies considerably from metal to 
metal, ranging generally from 0.01 mm. to 2 or 3 mm. and depends to a great 
extent on the previous mechanical and thermal treatment. The properties of the 
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metal aggregate, assuming it is free from coarse defects such as blow holes and 
cracks, depend on a combination of the properties of the individual crystals treated 
singly, the properties of the boundaries between neighbouring crystals (and of any 
alloying material precipitated there), and on any modification arising from the 
constraint of one crystal by its neighbours. 


A single crystal of a metal consists of atoms regularly arranged in a three- 
dimensional network. The distance between neighbouring atoms is of the order of 
two Angstrém units (=2 x 10-* cm.), and thus each crystal in a metal aggregate is 
many thousand atoms across. This assemblage of atoms acts as a three-dimensional 
diffraction grating for electro-magnetic waves, of wavelengths of the order of one 
Angstr6m unit. Such electro-magnetic radiations are called X-rays which may be 
produced by bombarding elements with electrons. The X-radiation produced by 
such a bombardment consists of a continuous range of wavelengths, i.e. “ white ” 
X-rays, with a superimposed line spectrum of X-rays produced by the quantum 
emission of radiation when an excited atom in the bombarded element falls to a 
lower energy state. In the X-ray diffraction work discussed here it is essential to 
know the wavelength of the radiation used and we are interested therefore only in 
the X-rays of this line spectrum. These are allowed to fall on the crystal under 
investigation and, since their wavelengths are known, measurements of the 
diffraction angles then enable the inter-atomic spacings of the crystal to be 
calculated accurately. 


It is impossible to deal in this paper with the complete theory of crystallographic 
notation and of X-ray diffraction in crystals, and reference should be made to a 
standard work’), A rough picture of the process can be gained from a consideration 
of the section through the assemblage of atoms shown in Fig. 1. Many series of 
imaginary planes, arranged in parallel groups, can be visualised as running through 
the three-dimensional atomic lattice; two such groups are shown (in both of these 
the planes are perpendicular to the paper). In each group the planes are equally 
spaced. To distinguish various planes a system of notation is used. The crystal is 
referred to three principal axes, parallel to the principal directions of symmetry of 
the crystal; many metals have a cubic symmetry and in this case the axes are 
orthogonal. Since the crystal structure is regular it can be built up by placing many 
small units side by side. The smallest unit that can be used has its sides parallel 
to the reference axes and of lengths a, b and c. A plane is designated (hkl) 
when the points at which it intercepts the axes are at distances a/h, b/k and c/l 
from the origin. If we wish to refer to the whole series of possible planes in the 
crystal parallel to (hkl), as we do in X-ray diffraction, the brackets are omitted, 


namely hkl. The planes shown in Fig. 1 would be 020 and 1 10. 


If the crystal is irradiated by a parallel beam of X-rays. of wavelength A, the 
X-rays scattered by each individual atom will, in general, be out of phase in all 
directions and there will be no diffracted beam. But if the incident X-rays make 
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Diffraction of X-rays by a crystal 


an angle @ with a particular set of planes spaced a distance d apart such that 
2d sin 6=A, then the X-rays scattered by each atom will be in phase in one particular 
direction and a diffracted beam will emerge from the crystal. The direction of this 
diffracted beam will be as if the incident beam had been reflected by the crystal 
planes. If the crystal is turned through all positions relative to the incident beam, 
diffracted beams will be produced by the various planes with spacings d,, d, . . . 
at angles 6,, 0, . .. such that the Bragg equation given above is satisfied in each 
case. The same will occur if, instead of turning one crystal to all positions, we 
irradiate many crystals randomly orientated in an aggregate. Thus if a narrow beam 
(of the order of 1 mm. cross-section) is directed at a polycrystalline metal, the 
diffracted X-rays will lie on the surfaces of a series of cones of semi-vertical angles 
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26,, 20, . .. , and with the incident beam as axis, each cone corresponding to 
diffraction by a particular set of Akl planes. The diffractions are usually recorded 
on X-ray films placed so as to intersect these cones; in the usual back reflection 
technique the film is a strip placed perpendicular to the incident beam and the 
diffractions are recorded as arcs of circles. 


If the crystals in the aggregate are larger than 10-* cm. and the atoms are 
regularly arranged so that the spacings in all parts of all crystals are the same, then 
the diffracted X-ray beams will be sharply defined. After certain treatments, 
notably “cold working,” broadened diffraction lines are obtained from the metal, 
and this can be attributed either to a small grain size“) or to variations in 
lattice spacing“). 


If the angles 6 corresponding to the maxima of the diffraction line are measured, 
then d, the average interplanar spacing of all the crystals contributing diffracted 
beams to that line can be determined. If d is changed by stressing the lattice, 
6 will change. This is the principle used in measuring strains; the change in @ 
is measured and the strain of those particular lattice planes is then obtained from the 
differential of the Bragg equation, namely 


= — cot 680. 


It is apparent from this equation that the nearer @ is to 90°, the larger is 66 for a 
given strain. In actual fact, the strains experimentally observed (which are rarely 
greater than the limits of elastic strain mechanically determined on large specimens) 
are so small that 6 must be greater than 75° if 56 is to be measured experimentally 
with sufficient accuracy. 


One point which cannot be emphasised too strongly regarding the diffraction of 
X-rays by a metal aggregate is that only certain crystals in the aggregate, those with 
the normals to the diffracting planes inclined at an angle 90°—@ to the incident 
beam contribute to the diffraction line. Measurements on a particular diffraction 
line give information about these particular crystals only. In Fig. 2, which illustrates 
this point, only the shaded crystals are reflecting. It should be noted that the 
crystal may have any orientation about the normal to the diffracting planes. This 
both limits the technique and makes it more valuable. 


Another point worth noticing is that X-rays are quickly absorbed by metals, 
and in general more than 50 per cent. of the diffracted intensity comes from a surface 
layer of less than 0.02 mm. in depth. Thus it is only the surface of the aggregate 
which is being examined. 
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METAL AGGREGATE 


NORMALS TO DIFFRACTING PLANES 
Fig. 2 
The selective action of the X-rays. Only certain grains in the aggregate diffract 


3. EXPERIMENTAL METHOD OF MEASURING STRAINS 


Extremely full and excellent accounts have already been published: * 7) of 
the type of apparatus used and the methods employed. These are all very similar, 
the difference between the methods employed by various authors being largely in 
the technique used to eliminate experimental errors in the determination of 6 values. 


The method employed by the author’) is indicated in Fig. 3. The wavelength 
of the X-rays used (the possible ones are, in general, limited to Kx emission lines of 
Cr, Mn, Fe, Co, Ni and Cu) is such that a diffraction line is obtained with 6 near 
80°. A narrow pencil of these rays falls on the specimen and the diffracted beams 
are recorded on a flat strip of film. On the surface of the specimen is placed a thin 
film of annealed silver powder, which also gives diffraction lines recorded on the 
film; these serve to calibrate the film-specimen distance and are used in the 
elimination of other experimental errors. The value of @ for the diffraction line from 
the specimen can be determined from the photograph, a typical example of which 
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COBALT TARGET 


\X-RAY FILM 


Fig. 3 
Schematic diagram of experimental arrangement 


is shown in Fig. 4. It should be noted that the Co Kz radiation used actually consists 
of two wavelengths 2, and «, very close together, and thus each diffracting plane 
gives two diffracted lines on the film. 


It is possible to examine the variation of strain with depth below the surface 
by etching off successive layers and taking an X-ray exposure from each freshly 
exposed surface. The depth of penetration of the X-rays is so small that fresh 
material can be examined after etching off a layer as thin as 0.1 mm. It is important 
to choose a satisfactory etchant which removes the metal evenly with little inter- 
granular attack, and which does not produce lattice spacing changes through causing 


SILVER. ALUMINIUM 


Fig. 4 
Typical X-ray film obtained with apparatus of Fig. 3 (actual size) 
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hydrogen to diffuse into the metal; electrolytic etching methods seem to be worth 
a trial, but the author has found that in some cases these will change the parameter 
of strain-free metal. 


It is usual in the examination of a given region of the specimen to take three 
exposures, one with the incident beam normal to the surface, and two with the 
incident beams making an dngle of 45°—60° with the surface. For the last two 
photographs the incident beams lie in mutually perpendicular planes each containing 
the normal to the surface and a direction of principal stress in the specimen. If 
the metal of the specimen is assumed to be elastically isotropic, and if it is assumed 
that all the strain is due to large scale stresses which are uniform within the volume 
of specimen irradiated by the X-rays (two assumptions which are commonly made) 
then the differences in strain measured in the three directions can be interpreted in 
terms of stresses by a formula.” *.7) An advantage of the three-exposure technique 
is that it eliminates the necessity for determining the strain-free lattice spacing. 


The current literature rarely investigates either assumption properly. In general 
neither is true, and in some cases the results obtained by using them can be very 
far from true. The trouble arises from the fact that metal crystals are not isotropic, 
i.e. the value of some properties depends on the direction in which they are measured. 
For example, in the case of copper, the value of Young’s Modulus measured 
perpendicular to the (100) plane is twice the value measured perpendicular to the 
(111) plane. If the random aggregate of a sufficiently large number of such crystals 
is examined the aggregate as a whole is found to be approximately isotropic, but if 
only grains of certain orientations are selected out of the aggregate and examined 
then the average behaviour of these will not be isotropic. The rules of elastic and 
plastic isotropy also fail for aggregates as a whole if such aggregates show preferred 
orientation, i.e. all the crystals tend to have the same orientation. Boas‘) quotes 
cases where Young’s Modulus and the yield strengths of commercial rolled sheets 
of steel and of copper differ greatly when measured perpendicular and parallel to 
the rolling direction. 


The remainder of this paper is devoted to the examination of the effect of the 
elastic and plastic anisotropy on the interpretation of the strains measured by X-rays 
in terms of stresses. In this it is necessary to distinguish between macroscopic 
stresses, which are large scale body stresses remaining uniform throughout regions 
at least as large as that irradiated by the X-rays, and microscopic stresses which 
are only uniform throughout a region of the order of size of a single grain and which 
may vary rapidly from grain to grain. 


4. THE EFFECT OF ELASTIC ANISOTROPY 


In general, the strain « measured by X-rays in a polycrystalline aggregate in 
any direction relative to a stress oc, is related to the stress by a constant k, such that 
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e=ko. If the strains are within the elastic limits this will always be true whether 
the volume under examination is behaving isotropically or otherwise. The general 
method of determining « involves measuring ¢, and ¢, in two directions in a plane 
containing the normal to the surface and the stress; then o=(e,—e,)/(k, —k,). 
(The origin of measurement is so close to the surface that no stress can exist 
perpendicular to the surface.) Therefore one must know the value of the constants 
k appropriate to the direction 1, 2, 3, ... 


These can be determined experimentally by applying stresses to the metal 
considered and measuring the resultant strains by X-ray methods, employing the 
same diffraction line as will be used in later stress determination work. Such 
experimental determinations: °°) show that k differs from that which would be 
expected from isotropic considerations and, moreover, the value for k in a given 
direction depends on which diffraction line is examined. The variations can be as 
large as 25 per cent. for steel. This is conclusive proof that the departure from 
isotropic elasticity is considerable and must be taken into account. 


Since it is a laborious task to determine the value of the constant k experi- 
mentally, investigations have been made into the possibility of calculating k from 
the known elastic constants for single crystals of the metals. The partial success of 
the calculations is of interest because it gives an insight into the behaviour of 
differently orientated contiguous crystals when a stress is applied. The methods 
were first applied by Reuss’) and had no relation to X-ray work; Neerfeld‘*) used 
the methods but did not publish the formule, and ihe author®?) has recorded a 
full discussion. 


The fundamental difficulty to be overcome in the calculation is to allow for the 
inter-action of one crystal with its neighbour. Suppose two contiguous crystals with 
different orientations are both stress and strain free. Since they are anisotropic, 
when an equal parallel stress is present in each the strains will not be the same, and 
if they could deform independently there would be a discontinu'ty of strain at the 
boundary; alternatively, if each is strained in the same manner \here would be a 
discontinuity of stress at the boundary. It is not possible to preserve the continuity 
of both stress and strain across the boundary during the deformatio.u. Mathematically 
the problem is tackled in two steps; first the value of k is calc.tated assuming that 
the strains in all the crystals are the same, and then a separate calculation is 
performed assuming the stress is the same in all crysta’s. It has been shown 
experimentally that the value of k determined by X-ray methods is approximately 
the mean of the two calculated values, and this is taken as being sufficient 
justification for using the method. 


It is thus possible to obtain, by experiment or by calculation, values for the 
elastic constants appropriate for the translation of the measured strains into large 
scale macroscopic stresses in cases where no plastic deformation has occurred. 
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5. THE EFFECT OF PLASTIC ANISOTROPY 


It is useful, before discussing the effect of the plastic anisotropy of crystals, to 
consider the behaviour to be expected if the crystals were isotropic from the point 
of view of plastic deformation. 


Suppose experimental conditions are arranged such that strains are measured 
in a direction perpendicular to the applied uni-axial stress, say a tension. Up to the 
elastic limit the strain will be proportional to the applied stress. If a decrease in 
the interplanar spacing is denoted as a negative strain, then the gradient of the 
curve of stress against lattice strain will be negative and equal to v/E where v and 
E are the values of Poisson’s ratio and Young’s modulus appropriate to the particular 
diffraction line investigated. This is shown in Fig. 5a; Y represents the limit of 
elasticity. When the stress is greater than that corresponding to Y, plastic extension 
occurs and the external dimensions contract much more per unit increase in stress 
than in the elastic range. But this plastic deformation takes place by glide”) 
which does not directly change the inter-atomic spacing; this is only affected in so far 
as the crystal becomes stronger and can support a higher stress. It is to be expected 
that the curve of stress against lattice strain would continue beyond Y along the 
same straight line to fracture. If, in actual fact, it did so the method of interpreting 
Strains in terms of macroscopic stresses by a formula, allowing for the change in 
elastic constants, would be justified. 


Experiments to investigate the curve of stress against lattice strain have been 
undertaken by many workers, but perhaps Smith and Wood) were the first to 
demonstrate the important point most clearly. It was shown that the form of the 
stress—lattice strain curve was similar to that shown in Fig. 5b; the lattice strain 
became almost constant when the stress was increased beyond the limit of elasticity, 
and if a plastically extended specimen was unloaded there was a residual lattice 
strain OA, OB, etc., remaining in the specimen. 


The explanation of this residual latice strain has perplexed workers in the 
X-ray strain measurement field for some time, and even now it is not completely 
understood. The fact that the grains examined have strained lattices indicates that 
they are stressed and, since the stresses across any section through the specimen 
must balance, this stress must be compensated by an opposite stress in some other 
part of the specimen. 


The author has shown’) that the residual lattice strain observed after a given 
plastic deformation depends on the diffraction line examined. Some diffraction 
lines show positive residual strains while others show negative values; it thus appears 
that after the plastic deformation of a polycrystalline metal aggregate a microscopic 
system of stresses is set up in which some grains are in tension and some in com- 
pression. It has been shown further that this arises because the yield tension of a 
crystal depends on its orientation, and a successful quantitative treatment has been 
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APPLIED APPLIED 
TENSILE FRACTURE , TENSILE 
STRESS STRESS 


FRACTURE 


0 OAB + 
LATTICE STRAIN LATTICE STRAIN 
(a) (b) 
Fig. 5 
Stress—lattice strain curves 


(a) Hypothetical, for plastically isotropic crystals 
(6) Experimental, crystals are plastically anisotropic 


based on the principle that soft grains are driven into compression by neighbouring 
hard grains after a plastic extension of the aggregate. Unfortunately, it appears 
ihat this is not the explanation of the whole of the residual lattice strain, and further 
work is necessary to elucidate the matter completely. 


The full significance of the existence of the residual lattice strain and its 
associated microscopic stresses on the interpretation of strains measured in a 
specimen by X-ray diffraction methods is not yet fully appreciated. It is clear, 
however, that the usual method of interpretation can only be employed to determine 
macroscopic stresses in regions where there has been no plastic deformation, and it 
is consequently limited to castings, quenched components, welded structures, and 
so on. In cases of doubt, the strains in the specimen should be determined twice 
using two different diffraction lines, which it is known from previous experiments 
show widely different residual lattice strains after plastic deformation. If the two 
sets of strain measured in the unknown specimen are similar, after allowing for 
elastic anisotropy, then they can reasonably be interpreted in terms of a macroscopic 
stress. If they are different, microscopic stresses varying from grain to grain exist 
and, at the moment, rigid interpretation is not possible. 
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The conclusion that micro-stresses exist in plastically deformed aggregates js 
not new; it has been reached before following observations on the Bauschinyer 
effect and other macroscopic phenomena*), But their quantitative estimation js 
new and this may lead to further interesting research results. 


6. CONCLUSIONS 


It is possible to measure strains in certain crystals whose orientations are fixed 
by the direction of the incident X-ray beam aad the diffraction line examined. The 
difficulties are entirely experimental. If the region examined has not been plastically 
deformed it is permissible to interpret these strains in terms of macroscopic stresses 
which remain approximately constant within the volume examined by using one of 
the usual formule embodying the appropriate elastic constants. If plastic 
deformation has occurred in the region examined, micro-stresses which vary from 
grain to grain probably exist in addition to macro-stresses and then it is not possible 
to do more than make a rough estimate of their relative magnitudes. 


Such conclusions are at least moderately satisfactory to the operator of the 
X-ray technique, but are unlikely to satisfy the engineer who is anxious to use the 
components examined. He is more interested in the practical interpretation of 
these conclusions made by the X-ray operator. This second step is even more 
difficult to take since the X-ray method itself is still only in its development stage, 
but it is an essential part of the work. 


At the moment, the only work done has been to relate macroscopic stresses to 
engineering properties. This has been successful in several cases, notable examples 
being given by Frommer and Lloyd‘). One of these was the use of the X-ray 
technique to examine quenching stresses in an engine crankcase; trouble experienced 
during the quenching of these was entirely eliminated after a change of design 
suggested as a result of the X-ray stress measurement. Again, considerable improve- 
ments have been made in fatigue properties by producing compressive macroscopic 
“locked up” stresses in the surfaces of components subjected to fluctuating load 
(although here the X-ray measurements are not entirely free from the suspicion 
that some microscopic stress was interpreted as a macroscopic one). In general, it 
seems to be advantageous to have a compressive stress on the surface of the 
component balanced by an internal tension. 


So far no use has been made of the measured microscopic stresses. The fact 
that two contiguous grains are in different states of stress indicates that there must 
be a high stress gradient at the boundary. This leads to a high strain energy there 
which probably has a marked effect on corrosion rates, particularly of the 
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inter-granular type. It is possible too that microscopic stresses also influence fatigue 
properties greatly, but much more research is required before the relation of these 
stresses to other metallic properties is understood. 


While the help that can be derived by the engineer from X-ray strain measure- 
ments is at the moment tentative and inconclusive, there is a reasonable hope that 
much may be accomplished in the future. 
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BUCKLING OF A LONGITUDINALLY 
STIFFENED FLAT PANEL 
by 


H. L. COX, M.A., F.R.Ae.S. 
and 
J. R. RIDDELL, B.Eng. 
(Engineering Division, National Physical Laboratory) 


SUMMARY 


(a) Purpose of Investigation. To determine the least size of stringers necessary 
to prevent overall buckling of a flat stiffened panel before buckling of the plates 
between stringers. 


(b) Range of Investigation. The conditions for buckling of a longitudinally 
stiffened flat panel are established on the assumption that the rotational stiffness of 
the stringers is negligible. The results are applied to determine the limiting 
characteristics of the stringers to ensure that these members remain straight up to 
the stress at which the plates buckle between stringers. The analysis has been 
carried through in detail for panels with one, two, or three stringers and is capable 
of extension to four, five or more stringers. This extension appears unnecessary, 
because the three stringer case is moderately close to the limiting case of a wide 
panel with an indefinitely large number of stringers. 


(c) Conclusions. The stiffening effect of the stringers, having each the area 
of section A, and modulus of section A,k*, when attached to a sheet of thickness ¢ 
and length a, at spacing b depends upon the three ratios (k/t), (a/b) and (A,/ bt). 
In general, the size of stringer necessary to ensure that the stringers remain straight 
until the plates buckle between stringers may be represented by the conditions 


A, r 
bt and (=) (where > 

For a panel stiffened by one stringer only suitable values are A= 4.5 and «= 12.5, 
independent of the value of (a/b), while v=0.366 (a/b)? or 23 if 8 > (a/b)>6. 


Paper received November 1948 
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For two or more stringers the same general conditions hold, but it is impossible to 
assign values to A and » which would be fully satisfactory for all values of the 
ratio a/b. In these cases, therefore, it is more satisfactory to select the values of 
A and » appropriate to each value of the ratio a/b. These values may be obtained 
from Figs. 2 and 3 where curves of bt/A, against (k/t)* for varying values of (a/b) 
are shown for two and three stringers respectively. For an indefinitely large number 
of stringers the corresponding condition is 


(a/by 


and the second condition (k/t)’ > 0.366 (a/b) is virtually included in the first. 


(d) Further Developments. In the present analysis the resistance of the stringers 
to rotation is neglected. Inclusion of this characteristic in the basic analysis is 
straightforward, but practical interpretation of the results becomes much more 
difficult. Moreover, the effect of this stiffness in modifying the condition for buckling 
between stringers is likely to affect the conclusions more than the direct effect of 
their torsional stiffness in modifying the condition for buckling of the stiffened panel. 
This consideration is particularly important in cases of light stringers and heavy 
sheet, when the effective rotational stiffness of the stringers at buckling of the plates 
between stringers may be greatly reduced as a result of the end loading and may 
even be negative. Representation of the effect of both longitudinal stiffeners 
(stringers) and lateral stiffeners (ribs), including their stiffnesses to rotation is also 
straightforward; but interpretation of the results in forms suitable for practical 
application would be extremely difficult. 


Notation 
a length of panel 
t thickness of sheet of panel 
b stringer spacing across panel 
A, area Of section of one stringer 
k radius of gyration of one stringer, including the effect of attached 
sheet 
E Young’s modulus 
oc Poisson’s ratio (taken as 0.3 in calculations) 
E’ E/(1-o?) 
M positive integer such that Mb=total width of panel 
r, Ss, m, m’ positive integers 
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e uniform compressive strain on panel 
t\? 

(5) 
c rMb/a 
X 3 (1 —o*) (k/ty 
Y A,/ bt 
6? M?e/e, 
[ce (c+26)]: in Appendix I (only) «=2r/a 
[c (c—28)] 
= * functions of « and 6 (Section 7) 
D, D, 

f(3, 1) functions of « and 6 (Section 9) 
R Mé/xcY (Xc? — 6) (Section 9) 
u,v, Ww amplitudes of buckle of three individual stringers (Section 9) 
U,V, W functions of u, v, w (Section 9) 
te stress in plate along length of stringers (Appendix I) 
fy stress in plate transversely across stringers (Appendix I) 
q shear stress (Appendix I) 
b’ effective value of b 
w > = = W,, sin = sin ws =deflection of plate 
x,y rectangular co-ordinates 
v lateral displacement (Appendix I) 
h depth of stringer 
P, Q,H, K numerical constants 
constants 
Wr Fourier coefficients 


1. INTRODUCTION 


In the design of stiffened panels to carry compression it is necessary to estimate 
the stress at which the panel will buckle as a whole. When the sheet cover is thin, 
and the stiffeners widely spaced, so that the plates between the stiffeners buckle at 
a low stress, the conditions for overall stability depend almost entirely on the 
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characteristics of the stiffeners. The condition for overall instability of a grid of 
ribs and stringers is known, ?) and this analysis may be applied to the case of a 
covered grid (that is, a panel), when the plates are already buckled between the 
stiffeners (Ref. 3, Part IV); but the method is scarcely adequate to discussion of 
cases in which the plates between stringers are either unbuckled or only just buckled 
when overall instability occurs (cf. Ref. 3 at the foot of page 29). 


The condition for instability of a panel stiffened only laterally is also known, *) 
and the method of analysis used in this reference is directly applicable to a panel 


stiffened only longitudinally or both longitudinally and laterally, when the sheet. 


between stiffeners remains unbuckled up to the stress at which overall 
instability occurs. 


On the other hand, the effect of longitudinal stiffeners is considerably more 
complicated than that of lateral ones, because their effective flexural and torsional 
stiffnesses are functions of the end loads which they carry. Despite this complication, 
complete solution of the general problem may be feasible; it is not attempted here, 
chiefly because the practical meaning of the result would be extremely difficult to 
elucidate. For example, equation (19) of Ref. 2 is already fairly complex; yet, in 
effect, the torsion constant J is also a function of the variable #, because the effective 
torsional stiffness of the stringers depends on the end load which they carry. Were 
the complete solution for the buckling of a flat panel stiffened both longitudinally 
and faterally available, it would certainly be necessary to simplify it drastically 
before conclusions of practical value could be drawn. Therefore, it appears 
reasonable to effect this simplification ab initio by restricting attention to the case 
of a panel stiffened longitudinally only, and by ignoring the torsional stiffness of 
these stringers. 


Extension of the analysis to the panel stiffened both longitudinally and laterally 
involves certain difficulties in respect of mathematical summations, and the exact 
solution may be inaccessible on that account; if only an approximate solution were 
available the difficulty of interpretation would be still greater. In respect of normal 
constructions with stringers spaced fairly closely and with ribs wide apart, the 
present analysis, or a slight extension of it, may be adequate; and it appears that 
stringers spaced much more closely than ribs generally represent an efficient 
structure (cf. Ref. 3). 


Extension of the present analysis to include the effect of torsional stiffness of 
the stringers appears more necessary, particularly in respect of constructions in 
which light stringers tend to become torsionally or locally unstable before the plates 
buckle between stringers. In this case no special difficulty is likely to be encountered 
in the mathematical analysis; but interpretation of the results may be 
extremely tricky. 


The problem of the buckling of a longitudinally stiffened plate has previously 
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been treated by Timoshenko’. The present treatment follows the same lines and 
indeed differs only in that it proceeds to a solution in closed form, whereas 
Timoshenko used a process of successive approximation. 


The advantage of the closed form of solution lies not so much in improved 
accuracy as in easier interpretation of the nature of the solution. Briefly, the results 
show that the stringer has to fulfil two requirements: first that it should itself be 
stable in flexure up to the required limit: and second that it should be stout enough 
to resist the tendency of the plate to drag it into a plate buckling mode. The first 
requirement for a given form of stringer section puts a limit on the ratio of the 
depth of the stringer web to the thickness of the plating: the second puts a limit 
on the stringer area of cross-section or in effect on the stringer wall thickness. 


2. STATEMENT OF PROBLEM 


The panel of thickness f, length a, and total width Mb is assumed to be divided 
into M plates each of width b by (M — 1) stringers, having each the area of section 
A, and effective modulus of section A,k*. In estimating k’ for stringers attached at 
one side of the panel only, the effect of the attached sheet is to be taken into 
account (see Section 5). 


It is assumed that all four edges of the complete panel are simply supported. 
It is assumed further that the stringer offers no resistance to rotational movement of 
the attached sheet about the base of the web of the stringer. Over long wavelengths 
the actual resistance of the stringer to this rotation is associated mainly with its 
torsional stiffness, including the torsion-bending effect; over short wavelengths this 
resistance is associated more with local deformation of the stringer section and it 
is also influenced markedly by shear lag effects. In constructions in which the 
stringer spacing is no more than two or three times the stringer web depth, the 
deformation of the stringers associated with buckling of the plates between stringers 
is likely to be mainly of the local buckling type, and in practice the effective stiffness 
in rotation may be positive but small, or negative and moderately large; the latter 
case will occur when the stringer wall thickness is small, so that the stringers tend to 
become locally unstable before the plates buckle between stringers. 


The assumption that the resistance of the stringers to rotation is identically 
zero provides probably a very fair approximation for all cases in which this stiffness 
is in fact positive, although it is not necessarily conservative; when the stiffness at 
the moment when the plates buckle between stringers is in fact negative, the present 
analysis provides only a rough guide. 


Relative to axes Ox and Oy parallel to the length and width respectively of the 
panel along one pair of edges, the deflection of the stiffened panel after buckling is 
represented by the general expression 


(1) 
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The whole panel is assumed to be subjected to a uniform compressive strain ¢ 


and the panel is assumed to be completely free to expand laterally). 


3. METHOD OF ANALYSIS 


That part of the total strain energy of the panel which om. on the 
deflection w is 


(1/2M E’* \{ (W.r+Wy,)? dxdy — 4Ete \| w,"dxdy + 


m=M-—1 
+ | ic dx—tEAg | ov.) dx | (2) 
ow _ 
where w, represents and w,.= Cte. 


(cf. Ref. 5 formule (1) and (2) where u.=e and v,= —ce). 


By insertion of the expression (1) for w this becomes 


(3) 


Dividing throughout by } EMabt and writing e, = / bY, this becomes 


(4) 
By differentiation with respect to W,,, we have 


b 


+a (3 W,, sin met) sin =0, (5) 
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1 (5 rMb 
{3 * rb) ~ eles 
(6) 


Writing rMb/a=c, A,/bt=Y and M?e/e,=6?, this 
formula becomes 


m —1 ae 

~ { W,, sin 

{ (s? +c’)? — 4c76? } 


Multiplying by sin = > S and summing over s, denoting = W,, sin a by W.., 
we have 
sin” sin ni sx 
m=M-—1 4 
m=1 s 4c26? 
sin (smz /M) sin (sm’x/M) 
The sum = (4+) — 426? 
ay m’)=/M }—cos{ s(m+m’)z/M } (9) 
(s? + (s? + 
where 2? + B?=2c? and «?8?=c" (c? — 46”), 
2? =c(c+26) and B?=c (c—28). 
COS S$, —COS Sb, & { COS 5¢,—COS 5H, COS SP, — COS SH, 
COS $$, — COS _ { cosh x (7 « (=~ $2) } (11) 
22 sinh 
\ 2 2 
= - (12) 
a sinh 
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provided that ¢, and 9, are both positive, that is provided that m { m’. Thus the 
sum (9) reduces to 


E { 78 (M—m)/M } sinh { ~Bm’/M } 


8cé 8 sinh «8 
sinh { (M—m)/M } sinh { zam’/M 
a sinh 


a is always real, but 8 may be imaginary; in that case, writing i8 for 8, the first 
term has sin for sinh everywhere but is otherwise unaltered. 


With the substitution of (13) for the s-series, equation (8) represents a set of 
(M — 1) equations, linear and homogeneous in W,». The condition for buckling is 
that the determinant formed by the coefficients of W,,, in these equations should 
be identically zero. Regarded as an equation in @ this condition is extremely 
complicated, because it is necessary to find the least value of 4 for any permissible 
value of c. The permissible values of c correspond to integral values of r in the 
formula c=rMb/a, so that each value of the ratio a/b needs to be treated separately. 


The present purpose, however, is restricted to determination of the limiting 
values of X and Y which afford no solution of the determinantal equation for 
6 > M (cf. Ref. 6) because when 6=M, the plates buckle between the stringers 
(the stiffness in rotation of the stringers being assumed to be zero). For this 
purpose, it is in the first place practicable and convenient to avoid reference to the 
ratio a/b by assuming c capable of continuous variation. This assumption would be 
strictly legitimate only if the ratio a/b were infinitely great, and this condition is 
debarred by other considerations. Yet the exposition gains in clarity by the approach 
adopted, and the effect of finite values of a/b can be considered later. 


4. THE SINGLE STRINGER CASE 


In the case of a single stringer, M=2 and m=1 only, equation (8) becomes 


and the limiting values of X and Y, when 6=M =2, are related by the formula 


1 FC (4_ { tanh { _ tanh { [ce (c+4)] (15) 


[ce (c—4)] (c+4)] 
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VALUES OF | 
Y As 5 
5 
” IN THIS RECION OVERALL 
INSTABILITY IS PREVENTED. 
OF 5O 100 200 250 
Valves of C VALUES OF x=s st)” 


Fig. 1 
Buckling of a longitudinally stiffened flat panel. Condition for stability of a single central 
stringer. Stringer area of section As: Modulus of section Ask*: Panel thickness ¢ and total 
width 2b 


When c>2+4/3, the last factor in (15) is positive, and no limit for 
Y is indicated unless X <4/c? <1; this is not a practical case. When 
2+ 73>c>2- ¥3, the first term of the last factor takes the tan form and, since 
in that range 2 > /[c(4—c)] > 1, this term is negative, and so also is the whole 
factor. Therefore a positive limit for Y is indicated only if Xc? > 4. By trial of 
values of c in this range, it is found that the lowest limits of 1/Y (that is, the highest 
limits for A,/bt) correspond to values of c in the range 2— /3(=0.268) up to 
about 0.380. These limits are shown in Fig. 1 as the envelope of the straight lines 
relating 1/Y and X for a succession of values of c from 2— 4/3 upwards; this 
envelope is terminated by the line corresponding to c=0.380 because the lines for 
higher values of c lie entirely above that for c=0.380. 


When c < 2— 3, the first term in the last factor of (15) retains the tan form 
but is now positive and, since this term is always greater than the second term in 
this factor, the whole factor is positive. Then no limit for Y is indicated unless 
Xc? < 4; but if Xc? < 4, the term Xc’ — 6 in (14) is zero for some value of @ between 
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0 and 2 and buckling must occur at this lower value of 6*. Thus in practice it is 
essential that X > a*/b?, because the least possible value of c=Mb/a (when r=!) 
=2b/a for the single stringer case. 


It is then unnecessary to consider the limitation on Y except at values of 
X > (a/b). Therefore, referring to Fig. 1, for values of a/b > 8 (say), for which X 
must be greater than 64, it is permissible to represent the limitation on -Y by the 
sloping line shown dotted; the equation to this line is Y =(12.2)/(X — 34). Moreover, 
for a value of a/b <6, the least permissible value of c=(2b/a) is > 0.333, and 
with this limitation the condition Y > (12.2)/(X — 34) again represents conservative 
design. In the intermediate range 6 < (a/b) < 8, the condition Y > (12.2)/(X — 34) 
may be slightly too liberal for values of X only slightly greater than (a/b). For 
instance, if a/b=7 and X =57, the stated condition requires Y > 0.53 whereas the 
true condition, obtained from Fig. 1 with X¥=57 and c=0.286, requires approxi- 
mately 1/Y < 1.6 or Y > 0.67. This laxity in the conditions may readily be covered 
by increasing slightly the required value of X in the intermediate range. This is 
illustrated in Table I where the limit X > 64 is applied over the whole range 
8 > (a/b) > 6. 


TABLE 1 
Effect of Finite Values of a/b 
a/b 12 10 8 2(2+¥/3) 7 6 5.25 5 4 


x> 144 100 64 64 (55.7) 64(49) 36 27.7 25 = 16 
r=1 0.167 0.2 0.25 0.268 0.286 §0.333 0.380 04 0.5 
r=2 0.333 0.4 0.5 0.536 0.572 0.667 0.76 08 1.0 


For the values of a/b > 10 the recommended limit on Y is never unduly 
conservative, and, since its greatest possible value over this range is 
(12.2)/(100 — 34)=0.18, it is unlikely to be very restrictive. In the lower range of 
values of a/b, however, it is sometimes too conservative; for example, if a/b=8, 
and X =64 it requires Y > 0.40 whereas, since the only relevant value of c is 0.5, 
the true limit from Fig. 1 is 1/¥Y <4.5 or Y >0.22 only. Yet over the range 
6 > (a/b)> 5 the recommended limit is moderately accurate. If, for a low value 


* In general, for all values of c, if the zero value of X°—@ should happen to correspond exactly 
with the infinity of the last factor of (14), which occurs when B= —1, that is, 2=4(c+1/c), 
buckling of the stringer as a strut would occur at exactly the same stress as that at which the 
plate would buckle in a single half-wave across the whole width 2b and down the whole 
length a. There is then no inter-action between plate and stringer, so that formula (15) is 
irrelevant; but see Section 5 as to the effect of another type of inter-action between stringer 
and sheet cover on the value of k 
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of a/b the condition Y > (12.2)/(X — 34) proves inconveniently restrictive, there is 
therefore no alternative but to refer to Fig. 1 in detail to determine the true limit 
of Y. Otherwise the conditions X > (a/b), that is, (k/t? > 0.366 (a/by 
or 23 for values of a/b between 6 and 8, and Y > (12.2)/(X —34), that is, 
A,/ bt > 4.5/ { (k/t—12.5 } represent conservative design. 


Table I helps also to emphasise the dual nature of this class of instability, 
which is more apparent at the larger values of a/2b, when the limit on Y (= A,/bt) 
results from a value of c corresponding to r=2 or higher. True flexural instability 
of the stringer is then covered by the condition X > (a/b), whereas the limit on 
A,/bt is to ensure that the stringer is sufficiently stout to prevent the panel forcing 
it to buckle in two or more half-waves. The same applies when r=1, but the 
process is less obvious when both panel and stringer are prone to buckle in 
one half-wave. 


It is worth noting that the type of buckling which requires a minimum value 
of A,/bt is analogous to wrinkling of sandwich panels. 


5. THE EFFECTIVE VALUE OF k 


The effect on the value of k of attachment of the stringer at one side of the 
sheet cover is important. The strain energy terms in formula (2) represent the 
flexural energy of the sheet cover and the diminution of the energy associated with 
the mean compression of the sheet by the bowing due to the buckle; but they do 
not include the energy associated with distortion of the sheet cover by bending of 
the stringer, except by appropriate choice of the value of k. 


By analysis of the stress system in a plate of width 2b subjected to shearing 
forces applied along its centre line, an expression may be found (Appendix I) for 
the effective width 2b’ of the plate according to the half wavelength (a/r) over 
which the shears are applied. The ratio of this effective width to b is 


b’ 4[1+ { sinh Qarb/a)/(2zrb/a) } 


+0) (arb/ay'+ (G+) cosh Qarb/a)+5—0° (16) 


If a—> 00, b’—> b, whereas if a—> 0, b’—> { 2/(34+¢)= } (a/r). 


6. EXAMPLE: ONE STRINGER 

For example, if a=6b, from Table I the appropriate value of c is 0.333, where 
r=1. So from equation (16), taking «=0.3, b’/b=0.907. For a Z-section stringer 
of depth A the value of k? is then about 


0.907 
0.907 +4 (A,/b1) } 


{0.15 +0.25 
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so that for a stringer one inch deep on a plate 0.08 inch thick 


(=) =23.44 = 


1+0.55 (A,/ 


Then the condition A,/bt > 4.5/ { (k/t)?—12.5 } demands A,/bt > 0.094; whereas 
a stringer } inch deep would need A,/bt >0.22 approximately. Thus, for a stringer 
attached down the centre of a plate 6 inches wide and 18 inches long the thickness 
of the stringer section would need to be at least 0.015 inch for the one inch deep 
stringer or 0.047 inch for the } inch stringer (taking one flange width of 
stringer = h/ 4). 


These computations assume perfect connection between sheet and stringer, 
so that where riveting is used some allowance for imperfect connection should 
be made. 


7. THE TWO STRINGER CASE 
By application of formule (8) and (13), the condition for buckling is 


30 
(Xc?— 6%) +fi=th., . ‘ ‘ (17) 
hes sin (278 /3) sinh (27/3) 
hh= 87142 cos 21 1+2 cosh 
(18) 
sin (78/3) sinh (7x /3) 


81142 cos 2 {142 cosh (272/3)}" 
or, writing = 3, the limiting values of X and Y are determined by the relation 
1/Y =(xc/9)(9—Xc*) Ff, + fr), 
with a=/[c(6+c)] and f=V[c(6-c)]* . (9) 


For stability at values of 6 <3, it is necessary that Xc? >9 for the least 
permissible value of c, which is 3b/a; that is X > (@/b)’, as in Section 4, and this 
condition applies generally, however many stringers are used. Therefore 9—Xc° 
is always negative and positive values of Y result only from negative values of 
f,+f.. As c increases from zero, (f,+f,) and (f,—f,) are at first both positive, and 
they increase uniformly to infinity, (f,+f,) at B=1 and (f,—f,) at 8=2. For values 
of 8 exceeding these limits up to 8=3, both are negative; but for 2<f <3, 
(f, +f.) is less numerically than (f,—f,). 8 cannot be greater than 3 if c is to 


* See note immediately after equation (13) 
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C=3-2/2=0-172 
0-4 


wn 


= 


ere) 200 300 400 500 
Fig. 2 
Buckling of a longitudinally stiffened flat panel. Condition for stability of two stringers 
equally spaced. Stringer area of section A,: Modulus of section Ask?: Panel thickness ¢ and 
total width 35 


\ 


remain real. The only limit is represented therefore by (f,+f.), which corresponds 
to the two stringers buckling in the same sense. 


The value of c corresponding to 8=1 is 3—2./2=0.172, and the corresponding 
least value of a/b (for r=1) is 3/(0.172)=17.5. So high a ratio of a to b may be 
unusual in practice, so that interest attaches rather to the higher values of c. 


sin (78/3) sinh 


(xB/3)—1) ~ Deosh (20) 


Writing Bi 


with 2= /[c(6+c)] and B= /[c(6—c)], by trial of values of c, the system. of 
straight lines relating 1/Y to X shown in Fig. 2 may be constructed. 


As for the single stringer (Section 4), it is again possible to state a general 
conservative limit for A,/bt(=1/Y), such as Y > 26/ { (k/t?—73} represented 
by the sloping dotted line in Fig. 2; and the reservation with respect to the X limit 
concerns only values of a/b in the range 17.5 to 18.2, for which we must have 
X > 330 in place of (a/b)’. On the other hand, values of a/b of 18 and over may 
be unusual, and in respect of values less than 12, the general limit may be far too 
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conservative. At the same time, since the values of A,/bt required are relatively 
larger than for the single stringer case, the condition may be less readily fulfilled, 
so that greater precision in its statement is desirable. For these reasons the detail 
of the diagram (Fig. 2) for the larger values of c and at lower values of X is probably 
of more practical value than the general conservative limit. The use of the diagram 
in this way is illustrated below. 


8. EXAMPLE: TWO STRINGERS 


For example, if a=6b, the least value of c=3b/a=0.5. The effective width 
of the sheet acting with the stringer in flexure of the stringer may be taken as 
0.905*. Then for a Z-section stringer of depth A the value of k? is about 


0.15 +0.25 h 
0.90+ 2(4,/b1) 


the factor 3 being in respect of two stringers and three plates. For a stringer one 


2 3 
inch deep on a plate 0.08 inch thick, (*) =23.44 While the 


line for c=0.5 in Fig. 2 demands (A,/bt)> 11.07/ { (k/t)’-—13.19}. Hence 
A,/ bt > 0.256; whereas if h=0.75, A,/bt > about 0.80 is required. 


It should be noted that, in this example, the general limit Y > 26/ { (k/t)? — 73 } 
is inapplicable, for (k/t)’ > 23.4+39=62.4, which leads to negative values of Y. 


9. THE THREE STRINGER CASE 
Denoting the function (13) by 3 f(M-—m, m’), and writing R for 


Mé/[xcY (Xc? —6)], we have for the three stringer case 


{R+f(3,1)} u+f (2, )v+f (1,1) w=0 
f(2,l)u+ {R+f(2,2)} v+f(2, )w=0 }, 
Du+f(2,)v+ {R+F(3,1)} w=0 | 


where u, v and w represent the amplitudes of buckle of the three individual stringers. 


* The case differs slightly from that worked in Appendix I; but no great precision in estimating 
the effective width is needed, so that the value previously computed in Section 6 is used 
again here 
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If we substitute 


ut+ /2v+w=4U u=U+V+W 
u— /2v+w=4V then v=/2(U-V) }. . (22) 
u-w =2W w=U+V-W 


Substituting (22) in (21) and using the relation f (3, 1)+f(1, 1I)=f (2, 2)t 
{R+f(2, 2)+ /2f(2, U+ 
+ {R+f(2,2)— 1)} V+ { R+f(, 1) } W=0 


[V2{R+f(2, 2)} +2f(2, 
+[2f (2, 1I)—- =0 


{R+f (2, 2)+ /2f(2, 1)} U+ 
+ D} V— } W=0 


Hence, either W=0 or R+f (3, 1)-f (1, 1)=0. 


(sinh 3Y—sinh Y)sinhY 
Now f (3, 1)—f (1, 1) has the form J sinh 4y a similar factor, 


and this first factor may be rewritten 


vsinh4y 2y 
which is identical with the form for the single stringer case (cf. equation (14)). 


This constitutes a complete solution of (23) if U=V=0; that is, if v=0, and 
u and w are equal and opposite; this is correct. 


The other two solutions with W=0 both make u=w; they are 
U=0 and R+f (2, 2)— /2f(2, 1)=0, 
or V=0 and R+f (2, 2)+ /2f (2, 1)=0. 


In the former case, the central stringer moves against the outer pair; in the 
latter all three stringers buckle in the same sense. It may be shown that the latter 
mode of buckling demands the highest standard of stringer stiffness. 


+ Because (sin 3¥+sin ¥) sin ¥=2 sin cos ¥ sin Y=sin? 
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Fig. 3 


Buckling of a longitudinally stiffened flat panel. Condition for stability of three stringers 
equally spaced. Stringer area of section A;: Modulus of section A;sk?: Panel thickness ¢ and 
total width 45 


The condition to be satisfied is therefore (for 6 < M=4) 
1/Y > (xc/16)(16— Xc’) { f (2, 2)+ /2f(2, 1)}, 


sin? (78/2) sinh? (r2/2) 
where FQ, 2)= Bsin(z8) 


1)= sin (78/2) sin (78/4) sinh (7/2) sinh (72/4) 


B sin (xf) x sinh (7x) 


a= /[c(8+c)] and B= V[c(8—c)]*. 


This condition is represented graphically in Fig. 3. Once again, it is possible 
to state a general limit for A,/bt(=1/Y), such as Y > 89/ { (k/t)?—243 } 
is represented by the sloping dotted line in this figure. However, this limit is 


, which 


* See note immediately after equation (13) 


ES 


BUCKLING OF A LONGITUDINALLY STIFFENED FLAT PANEL 


inapplicable for values of k/t less than 243 and, except for very large values of X 
with small values of c, it is more desirable to refer to the details in Fig. 3. 


10. EXAMPLE: THREE STRINGER CASE 


For example, if a/b=6, the least value of c=rMb/a= * and the condition 


for stability is Y > (33.9)/(X — 36) which is obtained by substitution in (24). For 
stringers one inch deep on a plate 0.08 inch thick we have 
= 23.44 the factor 0.83 being (0.90) for three stringers and 


1 + 0.83 (A,/bt)’ 
four plates with effective width 10 per cent. less than the full width. 


Then (A,/ bt) > (12.4)/ { (k/t) — 13.2 } 


> (12.4)/ { 10.2 + (39.0)/[1 + 0.83 (A,/bt))}, 


which demands A,/bt > 0.300; whereas stringers inch deep would require 


A,/bt > 1.07. 


ll. THE WIDE PANEL 


The solutions developed in Sections 3-10 and represented in Figs. 1 to 3 could 
be continued for four, five and more stringers; but it is not worth while, because 
the three stringer case is already moderately close to the limiting case of a wide 
panel with an indefinitely large number of stringers. In this limit, the effect of the 
restraint at the side edges of the panel may be ignored, and the requirement is 
merely to make the Euler load for the “ wide strut” (Ref. 3, p. 8) greater than the 
buckling load of the plates between stringers. The Euler load is 7*EA,k*/a* and 
the corresponding strain is =? (k/a)? A,/(A,+5t). The plate buckling strain is 
{ x*/[3 (1—o*)] } (t/b¥. Therefore the requirement is 


3 (1 —o*) (k/ay’ A,/(A,+bt) > (t/bY, 
XY/(Y¥+1)>(a@/b¥ where and Y=A,/bt. 


Writing this in the form Y > (a/b)’/ { X —(a/b)* } , the similarity to formule 
(15), (19) and (24) is apparent. 
Combining this equation with those obtained in Sections 6, 8 and 10 where 


241 


or 
| 


H. L. COX AND J. R. RIDDELL 


a/b=6 we have the general formula, applicable to any number of stringers, 
Y >»/ { where has the following values : — 


For one stringer 12.2 
For two stringers 30.2 
For three stringers 33.9 


For many stringers 36.0 


It is possible that » might rise above 36.0 intermediately. 
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APPENDIX I 


Effective width of a wide plate attached to a single stringer 


Taking axis Ox along the length of the stringer down the centre line of the 
plate and Oy transversely across one end and denoting the stress in the plate in 
direction Ox by f,, in direction Oy by f, and the shear stress by q, the general 
expressions for the stresses in either half of the plate corresponding to a harmonic 
variation of stress down the panel length are 


f-=(P cosh zy + Qzy sinh «y+H sinh zy + Kzy cosh zy) sin «x, 
q=— {(P—Q)sinh zy + Qzy cosh zy + (H — K) cosh zy + Kay sinh zy } cos «x, 


fy= — { (P—2Q) cosh «y+ Qzy sinh «y+ (H —2K) sinh «ay + cosh ay } sin «x. 


242 


| 
iz 
|| 


BUCKLING OF A LONGITUDINALLY STIFFENED FLAT PANEL 


— 4 _ 
b 
08 7 Z 
/ 
0-6 yy, 
0-4 
0-2 
0 5 10 @ 15 20 
rb 


Fig. 4 
Buckling of a longitudinally stiffened flat panel. Effective width of a wide plate attached 
to a single stringer 


From the relation Edv/dy=f,—of., we find the lateral displacement 


v= —(1/Ex)[ { l+c) } sinh Qzy cosh «y+ 
+ } cosh Kzy sinh zy] sin ax. 
Then since v=0 at y=0, (l+o)H=(3+0)K. . (Al 


The edges of the plate at y=b being free, we have, writing C for cosh«b, S$ for 
sinh 2b and B for «ab, and using (A1) 


. . (A2) 
—fy=PC+Q (BS—2C)+ K [ { (l—o)/(1+¢) } $+ BC]=0. (A3) 


Hence P and Q may be expressed also in terms of K. 
The average stress across the width of the plate is 


f. dy=(1/2) { . (A4) 
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B 
so that the effective width, b’, which is the ratio of | f. dy to the value of f, at 
x=0 is given by : 

(o’/b)= { (P—Q)S+QBC+(H—K)(C—1)+KBS } /PB. 
After substitution for P, Q and H in terms of K this reduces to 


4[1+ {sinh (2rb/a)/(2nrb/a) } 
~ cosh (2arb/a)+5-—0’ 


5’ 
b 


where «=zr/a has been substituted. 


The variation of the ratio b’/b with a/rb is shown in Fig. 4. 
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APPROXIMATE CALCULATION OF 
THE LAMINAR BOUNDARY LAYER 
by 


B. THWAITES, Ph.D. 
(Aeronautics Department, Imperial College, London) 


INTRODUCTION 


The steady two-dimensional flow of viscous incompressible fluid in the boundary 
layer along a solid boundary, which is governed by Prandtl’s approximation to the 
full equations of motion, presents a problem which in general is as intractable as 
any in applied mathematics. The problem, however, has such an immediate and 
necessary application that approximate methods of varying accuracy which go 
beyond the formal processes of expansions in series and so on, have been devised 
for the rapid calculation of the principal characteristics of the laminar boundary- 
layer, the variation of pressure along the surface being known. Such methods usually 
represent approximately the boundary-layer velocity distribution at any point by one 
of a known family of distributions whose spacing along the surface is determined by 
some means, often by the use of Kérmdén’s momentum equation. 


In the first main part of this paper, Sections 3-5, all known velocity distributions 
from exact and approximate solutions are collected and compared in a manner 
which shows clearly their potentialities or limitations when used as the bases of 
approximate methods. This critical comparison explains that it is possible for an 
approximate method to yield the exact values of the more important characteristics 
of a flow, and enables a method of calculation to be constructed in the second part 
of the paper, Sections 6 and 7. 


The momentum equation, of which a full discussion is given in Sections 1 and 2, 
is the basis of the construction, but is not necessary for the final working, of the 
new method. The method, which is shown to give good results and to be simple 
and speedy in application, can be used with confidence in regions of increasing 
pressure and, in particular, predicts separation of flow with good accuracy. 


Paper read before the VIIth International Congress of Applied Mechanics, Imperial College, 
London, September 1948 
[The Aeronautical Quarterly, Vol. 1, November 1949] 
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Notation 
Note: 


x,y 
u, v 


Vv 
U., p.. k 
A, A, 

f 


fA 
T(A), JA) 
l 
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U’, U” are derivatives of U with respect to x 
f f’, f” are derivatives of f with respect to 7 
F’, F”, F’’ are derivatives of F with respect to 7. 


distances parallel and perpendicular to the plane boundary y=0) 
velocity components parallel and perpendicular to the boundary 
the stream velocity 

the value of v at y=0 

the stream function 


oo 


(1 dy, the displacement thickness 


2 


u u 
“(1 dy, the momentum thickness 


a length of the order of 5* 
a chord length 

a parameter, variously defined, denoting distance perpendicular to ( 
the boundary 
6/8 
5* 

6 

the kinematic coefficient of viscosity 

various constants occurring in definitions for U 
numerical parameters 

a function of » defined in Section 4 (ii) 


certain functions of A 
U oy y=0 


6? (07u 
U 


an index of x/c 
2k 

k+1 

(k+1)) 

(K+ FQ) 

a suction parameter 

a function of Kn 


l 
| 
8 
c 
7 
0 
H 
| 
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L (m) 2[ (m)+2)m+l(m)] 
M a function of 6? and U 

an x-interval 

various coefficients 


F a function of various parameters according to context 
t 
U 
Ué 
Ro 
v 
Bix 
Bo 
certain sums 


1. PRELIMINARY DISCUSSION OF EXISTING METHODS 


(i) Prandtl’s equation of motion of steady two-dimensional boundary-layer 
flow is 
(u, v) being velocity components, (x, y) Cartesian co-ordinates, U the stream velocity 
just outside the boundary layer and y=0 the solid boundary. v is the kinematic 


viscosity. 


If equation (1) is integrated from y=0 to y=00, the first integral or momentum 
equation is obtained asf 


dé fou 


in which the momentum thickness @= j 7 (1 ~ x) dy, the displacement thickness 


(1 x) dy, H=6*/@ and the equation of continuity 0 is used. 


The equation of motion (1) becomes, when y=0 


a 


Of the methods which use equation (2) directly, Pohlhausen’s"? is the simplest, 
y which assumes a family of velocity distributions analytically defined. It is well 


+ It should be noted that, in the case of a porous surface through which there is a normal 


velocity Yo, an extra term z should be added to the right hand side of the equation 
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known to give poor results in regions of rising pressure. Various other authors have 
tried to improve and extend his method by assuming different families of 
distributions. 


It is possible to devise other methods for giving a correspondence between a 
known family of velocity distributions and a certain flow. For example, Howarth?) 
uses the distributions which occur in the solution of U =8, — 8,x as the basis of an 
approximate method for regions of rising pressure. Falkner‘) uses the distributions 
corresponding to the flows U =U, (x/c)* as a basis for an elaborate method in which 
the momentum equation is incidentally used. In Ref. 4, he simplifies this work and 
also suggests a form for @ in terms of U directly. Young and Winterbottom? also 
gave 6 as a function of the stream velocity U in a form which is close to that derived 
in this paper. 


The Karman momentum equation, being an exact integral of the equation of 
motion, should be able to yield good results but, without knowledge of any exact 
solutions, such an equation would be practically useless. It is one of the intentions 
of the present paper to show how known exact solutions can themselves be given 
accurately by correct use of the momentum equation. A method of solution is 
finally derived which is virtually the best possible under the general conditions of 
the use of the momentum equation. 


(ii) The assumption is usually made that to a sufficient degree of accuracy the 


velocity distribution — 7 (x, y) can be represented by the single infinity of 


distributions 


u y 


8 is some length of the order of the boundary layer thickness, A is a function of x 
and F is a function of two variables which satisfies the usual boundary conditions: 


F (0, 4)=0; F(z. a) = 1, y>Ab. 


If A is finite, the approximation assumes that the stream velocity U is actually 
reached at a finite distance from the boundary. 


Equation (3) imposes another condition upon F which is 


, 


2 
F” (0,A) or — =F” (0, A), (5) 


in which dashes on F denote differentiation with respect to y/5. (5) is an identity 
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which must be satisfied at all points of the boundary and which suggests that A 
can be conveniently defined by A= —8*U’/v. The significance of this relation is that 
it relates the thickness of the boundary layer and the stream-velocity gradient to A 
or to the “shape” of the velocity distribution. 


From (4), expressions are obtainable for 6, H and (*) in the form: 
y=0 


6 500), = H=H (A), in which 9, H, f are functions of A 


depending on the form of F originally chosen. By substituting these expressions 


in (2) and with A= —- — we obtain, after some re-arrangement: 
d (=) Of 
yr) 


which can be reduced to 


dx U” U’ \ 
dx = IW +J (A) U . . (6) 


1 (A), J (A) being determinable functions. 


(6) can be integrated to obtain the distribution of A with respect to x, whence 
all the other characteristics follow at once. 


(iii) The method given above is the standard method of using the momentum 
equation. In (6), A is the dependent variable—some writers have made 6? the 
dependent variable, which is a slight alteration but has no effect on the comments 
made below. 


First, to carry out the whole process from the beginning is tedious. A large 
amount of algebra is involved in obtaining © (A), H (A), f (A), 7 (A) and J (A). When 
(6) has been integrated numerically, further work is needed to find the distributions 
of 6, H and (“) with x. 


y=0 
In equation (6), the quantity U” cannot be evaluated with much accuracy from 
an experimental set of readings of the pressure distribution along a surface. Further 
there is a singularity in the equation at U =0 and it is necessary to find the correct 
value of A to be taken at a stagnation point. 


The introduction of 5 as a convenient but unrequired parameter involves not 
only general complication in the handling of the equations, but additional 
computation to obtain @ and H. 


. 

| 
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Lastly, there is always the possibility that the solution requires values of A for 
which the method is unworkable (for example, when A> 12 in Pohlhausen’s 
method), and also that for some values of 4 equation (4) no longer gives a reasonable 
approximation since F may have a maximum greater than unity. 


Next to be considered is the general approach to the equations (2) and (3). 
These equations are regarded as simultaneous equations in 6 and A, which are the 
parameters in terms of which the required variables are expressed. It happens 
that the second equation allows the very simple relation, A= — 5*U’/v between 6 and 
A and it is natural to make analytic use of this to turn the momentum equation 
into an equation for either 5 or A. 


The present author considers this to be rather the wrong approach to the 
equations in that it obscures the initial assumptions of a certain family of 
distributions and confuses the true relationship between the two equations (2) and 
(3). In what follows, a different approach will be explained and a simplified 
treatment of the momentum equation given. 


2. THE NEW APPROACH 


(i) In general, the intention is to avoid the introduction of unnecessary para- 
meters such as 3, to avoid all algebraic work (except, of course, in the actual! 
integration of the equation), to simplify to a great extent this integration, to ensure 
that the method can always be worked, and to avoid the necessity for knowledge 
of U”. These aims will be attained largely by the use of @ as the principal 
dependent variable. 


The main purpose of any approximate method is to give as accurately as 
possible the distribution of skin friction (and therewith the point of separation, if it 
exists) and of the momentum or displacement thickness. A good representation of 
the velocity distributions is not necessary. In the method to be presented, a 
synthesis of all known accurate solutions enables the first purpose to be achieved 
with the greatest degree of accuracy possible, and a good representation of the 
velocity distribution can also be obtained simply. 


(ii) The essential features of equations (2) and (3) are the two terms 


2u) 


2 
From the second equation v (3) is equal to —UU’. To integrate the 
y=0 


momentum equation the value of () is required. Thus the fundamental 
y=0 
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requirement is a relation between the two terms in (7). For if some such relation 


is assumed (e.g. by the use of a family of velocity distributions such as (4)) then 
( ) is known as a function of x since the value, — — , of (53) is known. 
\OV/y=0 y=0 

The two terms (7) indicate the behaviour of the velocity distribution at the 
boundary. Its general “shape” elsewhere is indicated by the value of H which also 
occurs in the momentum equation. This equation therefore is not only concerned 
with boundary values, but relates these with the velocity distribution as a whole. 
For this reason, the momentum equation is capable of giving good results. Thus 


0? 
a relation between H and (5*) is also required. 
y=0 


The construction and use of these required relationships is now considered. 


3. THE COMPARISON OF BOUNDARY LAYER VELOCITY 
DISTRIBUTIONS 


All known velocity distributions are now compared by classifying in a particular 
way the values of their first and second derivatives at y=0 which we have seen, 
following (7), are of particular importance. 


Ou _U (3) 


these forms being chosen so that the numbers / and m depend only upon the “shape” 
of the distribution of velocity, and not upon the thickness of the boundary layer 
or the magnitude of the stream velocity. We are to investigate the relation between 
| and m and to regard / and also H as functions, /(m) and H (m), of m. I[(m) and 
H (m) are therefore computed for all known solutions, exact and approximate, of 
laminar boundary layer flow which are enumerated and examined in the 
next Section. 


4. EXAMINATION OF KNOWN SOLUTIONS 


Brief remarks are now made on all the known solutions of the boundary layer 
equation. Details of the computations involved in obtaining /, m and H will mostly 
be omitted, but the results are given in graphical form. Pohlhausen’s velocity 
distributions and the corresponding method of solution are discussed in some detail 
and fresh light is thrown upon the well-known failure of his method in certain cases. 
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(i) Pohlhausen’s profiles are given by 


It is easy to show that 


UX ( x 5A? 
=— and 3 - ai): 


6? A ( 5A? \? 


from which 


m therefore has maxima or minima at 
A= — 37.79, —17.76, +12, +28.19, 
at which m takes the values 


respectively. 


Now(5*) (12+A), and therefore when A= -—12, separation of flow 
y=0 


occurs, and for A < — 12, >) <0. Thus we are only interested in the range 
y=0 
of A,A=>-—12. As A increases from — 12, m decreases from 0.157 until a minimum 
in m is reached atA= +12. (At this point, it is well known that the solution becomes 
unworkable.t) Values of m less than —0.095 are therefore unattainable except for 
much larger values of A. These larger values of A however cannot be used since 


in a solution they would involve a discontinuity in A from A= 12 to A= 34.79 (where 


again m= -—0.095) and therefore discontinuities in 6 and (3) » which are 
y=0 


{+ For example, in Modern Developments in Fluid Dynamics (edited by S. Goldstein), 
Vol. 1, p. 161 
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inadmissible. Furthermore, for values of A > 12, the velocity distribution given by 
(9) has a maximum greater than unity which in itself would make a solution in this 
region inaccurate. The range of A available is therefore —12 <A < 12, for which 
0.157 > m>-—0.095. 


It may be pointed out here that in the problem of the flat plate in a uniform 
stream, with a constant normal velocity through the plate, the velocity distribution 
tends, as x —> 00, to the asymptotic suction distribution, for which m= -—0.25f. 
It is clear therefore that the Pohlhausen distributions are very far from reaching 
this value of m and could not be expected to give satisfactory results in this problem. 
(See, for example, Ref. 14, Section 1.) 


The functions /(m) and H (m) are easily calculated and they are shown in 
Figs. 1-4. 


(ii) U=U, (x/c¥ 

Falkner and Skan“) have obtained exact numerical solutions of this flow which 
Hartree) has amplified on the differential analyser. For any value of k, it is 
well known that the distributions have a constant “ shape.” The equation of motion 
is transformed by writing 


u==-, p=(Uvx}f(n) 


Ox vx 
— 
(12) becomes 


F’’ +FF” =B8 ((F’¥ 
with the boundary conditions 


F(0)=0, F’(O)=0, F’(Y)—1 as Y— 


+ The asymptotic suction distribution can be written as “ =1-e-»/%, Therefore 


U 
(oe U whence m=-—0.25. See the footnote to Section 4 (v) 


yao 40 


253 


) 

= 

| 
| 

|_| 


B. THWAITES 


The upper limit for m for />0 is 0.068 and in this case k= —0.0904 and 
1 (n)=0, i.e. the stream velocity decreases in such a way that a separation distribution 
is just maintained. 


It is difficult to determine whether there is a lower limit for m. Calculations 
have not given values of m less than -0.10, since Hartree. pointed out that the 
practical value of such distributions is not very great. This is true for flow along 
solid boundaries but in view of the flat plate constant-suction problem for which 
a value of m= - 0.25 is required, it is clear that a greater range of solutions of this 
flow would be useful. 


(iii) Linearly decreasing pressure gradient: (U=8, — B,x). 


This flow has been treated by Howarth’) using expansion in series, the accuracy 
of which is fully discussed in the original report. Since U’ < 0, only positive values 
of m occur. Separation takes place when 8,x/8,=0.120 and m=0.084. 


(iv) Schubauer’s ellipse: Hartree’s solution"? 


Hartree has solved the boundary layer equations on his differential analyser for 
the pressure distribution which Schubauer observed on an ellipse. Near the forward 
stagnation point a series solution was applied and the analyser solution was started 
at a distance 0.067c from the stagnation point. Hartree actually did not find 
separation occurring before, or at, the observed point of separation and points out, 
in a full discussion, that it is meaningless to pursue the solution, using the observed 
pressure distribution, past the observed point of separation. He also points out 
how sensitive the solution is near the observed separation point, extremely small 
changes in pressure gradient being sufficient for the solution to give a separation. 
Howarth has discussed this experiment in relation to approximate solutions in Ref. 2. 


For Hartree’s solution, m takes the value 0.0784 at the limit of calculation. In 
the region U’>0, m takes the value - 0.0854 at the stagnation point which 
corresponds to the value k=1 in (12) and increases steadily to zero at the pressure 
minimum. 


This is the special case k=0 of equation (12). The boundary conditions are 


fO)=K, f()=0, f(~o)=1 


and the solution corresponds to the uniform flow past a flat plate through which 
x 


) . K>0 therefore corresponds to suction 


there is a normal velocity v.= — *( 
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and K=O gives the well-known Blasius’ distribution for uniform flow past a solid 


; boundary. The equation has been solved for various values of K—Ref. 9 deals 
with positive values of K and Schlichting”) gives solutions for positive and negative 
values of K. This family of velocity distributions is interesting for several reasons, 

S especially for K >0. Watson‘'") has proved that as K —> +00 the velocity distribu- 

> tion tends to the asymptotic suction distribution.t Thus the negative range for m 

3 is 0 => m >-—0.25 which covers the range necessary in the flat plate constant-suction 

problem. For K <0, the distributions are curious (physically this corresponds to 


) blowing out of the surface with velocity proportional to x~-!). As K increases 
m increases but reaches a maximum of about 0.031 and then appears to tend to 
zero as K tends to —°0. /(m) steadily decreases as K decreases and also tends 
to zero. 


(vi) Iglisch’s solution of the flat plate, constant-suction problem.“?) 


This and the next solution described are of considerable interest because they 
deal with velocity distributions obtained when boundary suction exists. There is 
reason to suppose that the general type of distributions under conditions of boundary 
suction may show rather different characteristics from those at solid boundaries, 
especially in regions where suction is delaying or preventing separation. This 
solution of Iglisch was obtained by numerical methods, there being no exact solution 
in finite terms known. The distribution is Blasius’ at the leading-edge of the plate“*) 
and changes progressively until the asymptotic suction distribution is reached. Thus 
0<m<-0.25. 


wh 20: 


This represents the flow from a stagnation point under conditions of boundary 
suction. The solution is a particular case of Falkner’s general solution (equation (12)) 
i with kK=1 and f(0)=K, and has been calculated for various values of K by 
Schlichting and Bussmann®). In this flow, the velocity gradient reinforces the 
general effect of the suction, and m takes only negative values. 


f(y)=K+ (Kn), f=0', f’=Ko”, f” =K*¢’”, equation (13) becomes +49” = oo”, 


with 9’ (0)=9 (0)=0, 9’ As the equation tends to 9” +4¢”=0, of which the 


solution is (t)=1—e-* or 


If U is not constant a similar result holds, for if U=o0 (vo), u=o(Vo) the continuity 
2 
equation gives v=v»)+O(U), so that the equation of motion is vo a -v cs =o (=) the 
solution of which as is =1—e" 
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5. COMMENTS AND COMPARISONS 


(i) In Section 4, all the known solutions of laminar boundary layer flow have 
been mentioned. The functions /(m) and H(m) have been drawn for all these 
solutions and are shown in Figs. 1-4. Each point on these curves represents a certain 
velocity distribution and therefore in these figures are collected a large range of 
different distributions, each of which is known to exist. 


Values of m divide naturally into the two ranges of positive and negative values, 
since for solid boundaries positive values of m occur when U’ <0, negative values 
when U’>0.+ (For conditions of continuous suction this is no longer true.) 


Consider first negative values of m. Fig. 1 shows the function /(m) drawn 
for the various flows. The range of m shown is 0 >m >-—0.12, which is sufficient 
for all distributions computed except for Iglisch’s solution. For this, as has been 
already pointed out, there is the point for which /=0.5, m= —0.25. It is remarkable 
that the various values of /(m) for any particular m do not differ by more than 
about 10 per cent., although the distributions are those of widely differing types of 
flow. The suction distributions for U=U, (x/c) differ by about 10 per cent. from 
the other distributions for — 0.02 > m >-—0.04, in which region these distributions 
differ from each other by only about three per cent. The general comment upon 


Fig. 1 therefore is one of surprise that the various curves of /(m) lie so closely 
together. Fig. 3 gives H (m) and here again the values of H (m) do not differ by 
more than five per cent. H (m) increases steadily with m, and its minimum value 
appears to be 2.0 at m= - 0.25. 


Figure 2 gives /(m) for m=>0, and it is at once obvious that no such close 
agreement exists as for m <0. We can dismiss the curve for f” +4ff” =0, for not 
only does this curve correspond to distributions obtained by blowing out through 
the boundary but also the shape of the curve itself has no practical interest. The 
three exact solutions are grouped fairly well together, with the solution for 
U=8,- 8,x in the middle. The value of m for separation in U =U, (x/c¥ is 0.068 
whereas the value for U=8,—,x is 0.084. Had Hartree predicted separation for 
Schubauer’s ellipse, the value of m would have been near to 0.084. The difference 
between the Schubauer and the Howarth curves lies in the fact that in the former 
a boundary layer of some thickness already exists where U’ becomes negative, while 
for the latter the boundary layer has zero thickness at m=0. Finally the Pohlhausen 
curve lies well away from the other curves. Pohlhausen’s /(m) takes values which 
are too large and this means that in regions of negative velocity gradient the skin 
friction is given values systematically too high. The value of m at separation is 0.157. 


+ For, from equation (3), 0=UU’+ > 
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The reason for Pohlhausen’s method giving inaccurate results in regions of 
increasing pressure is therefore clear, and it would seem that Falkner’s method, 
which is based on the profiles of U =U, (x/c)* would give better results. Fig. 4 
shows H (m) for m>0. Here also the various curves are fairly widely scattered, 
although those of the three exact solutions are roughly of the same general shape. 
Pohlhausen’s H (m) gives values which are much too low. 


(ii) In general, to any flow along a surface corresponds an /(m) curve and an 
H (m) curve. There is an exception in the case of flows with “similar” velocity 
distributions, when single values of /, m and H describe a flow. Thus in Figs. 1-4 
the curves for U =U, (x/c) do not represent a single flow but a very large variety 
of different flows and it is all the more remarkable that the curves lie closely to 
those of other flows. 


It was shown in Section 2 (ii) that integration of the momentum equation 
requires only a relationship between the terms in (7) or the /(m) function and also 
the H (m) function. Thus the probable accuracy obtained by using any approximate 
method upon a certain flow can be assessed by a comparison of the approximate 
and exact /(m) and H(m). Thus if an approximate method has the same / (m) and 
H (m) as an exact solution, then it will yield exact results. An immediate corollary 
is that, since no two exact solutions have identical /(m) and H (m), an approximate 
method which gives excellent results in one case may not be successful in another. 
In other words, there is a limit to the degree of accuracy obtainable by an approxi- 
mate method used on any flow. But, in fact, the accuracy obtainable in any case 
can be quite sufficient for practical purposes, and the formule developed in this 
paper give what are in effect the best results possible in all cases. 


6. THE CONSTRUCTION OF A METHOD 


(i) It was pointed out in Section 2 (ii) that the integration of (2) in conjunction 
with the boundary conditions (3) requires only the functions /(m) and H (m), and 
to construct a method of integration the choice of these functions has to be made. 


We refer now to Figs. 1-4, in which /(m) and H (m) are plotted for all known 
solutions. The choice of /(m) and H(m) for use in the present method will be 
made by taking an average value of these solutions. 


When m <0, the choice is easy, for the curves of Fig. 1 lie closely together. 
When m> 0 the choice is more difficult. It is first necessary to choose the value 
for m at separation. (The fact that an arbitrary choice is involved is the fundamental 
disadvantage of the use of the momentum equation.) For U=8,—8,x, m=0.084 
at separation and for U =U, (x/c)-°°***, 1=0 and m=0.068. For Hartree’s solution 
of Schubauer’s ellipse the value of m at separation, had it been predicted, would 
appear to be about 0.082. This last value will be taken in the present method to 
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represent separation, since the flow about the ellipse is similar to flow about an 
aerofoil to which approximate methods of calculation are most applicable. 


It is desirable to make /(m) and H (m) take certain values which correspond 
to particularly common exact solutions. Blasius’ distribution occurs at the leading 
edge of a plate in a stream if there is no stagnation point, and if there is a stagnation 
point the distribution there is that for U =U, (x/c). These two distributions should 
be included in our solution. The asymptotic suction distribution should also be 
included, if only because its large negative value of m will avoid the possibility of 
a breakdown in the method. 


Thus we require: 


U=U,: 1(0)=0.220 
H (0)=2.591 
Asymptotic suction distribution : 1(-0.25)=0.5 (14) 
H (-0.25)=2.0 
Separation distribution: 1 (0.082)=0 
H (0.082)=3.7 


U=U,x/c: 1(-0.0854)=0:359 
H (-0.0854)=2.218 


Furthermore the function /(m) appears to have an infinite gradient at /=0. 


It would be simple to express /(m) by a short analytic form. It is surprising 
how definite a curve one wants to draw when taking a reasonable mean of the 
curves in Fig. 1, but in fact it is difficult to choose a function for /(m) which is 
short, simple and suitable. Therefore we have no hesitation in defining numerically 
1 (m) for the present method, and in Table I /(m) is tabulated. This / (m) includes 
the points (14) but its value at m= — 0.0854 (corresponding to a stagnation point) 
1s 0.344 instead of 0.359, which is not a serious inaccuracy. 


The function H (m) is found in precisely the same way and is given in Table I. 
Values of / and H at values of m intermediate to those given can be found by 
simple interpolation. 


With /(m) and H (m) now fixed, the momentum equation can be integrated. 


(ii) The equations (2) and (3) become, after substituting from (8): — 


o=uu'+ 
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On re-arrangement and substitution the first of these becomes 


5 (H (m)+2)m+I1(m) 

and if L(m)=2[(A (m)+2)m+l(m)] ‘ GF 
2 192 

then =L(m), m= ; . & 


The method of solving equations (16) is now obvious. At any point of the 
system under consideration we know U, U’ and @. These define m and thus L (m) 
and therefore the gradient of 6 is known. A simple step-by-step method can be 
used and the only requirement for finding the distribution of @ with x is the function 
L(m) which is tabulated in Table I. For obtaining the other properties of the 
boundary layer, we require to know /(m) and H (m). 


It is clear that the aim to use @ as the principal dependent variable, to avoid 
algebraic work and to simplify greatly the equation to be integrated, has been 
attained. 


(iii) A method is required for the integration of (16) which is simple but 
which gives a high degree of accuracy. 


The equation is of the form 


2 


in which 6? and U are regarded as independent of one another. Let suffixes 0 and 1 
denote values at x=x,, x=x, respectively, and let x,—x,=5. We assume that 
x, U and partial derivatives of M of any order are O (1), and =o (1). 


dé? 
Then we have 6,7=6,? + +O (8) 
or from (17), 0,7=98,? + O (8) 
where + 5M (6,7, U.) : ‘ ; 


We wish to refine this formula to give @,? to O (8°). 
Now, 


M (6,7, U,)=M (8,7, U,)+O U,)+0 (8) (19) 
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In general, fr=fo+8 +0 (8) 


and +8 O (8?) 


and elimination of (f”), from these two equations gives 


i 
fi=fot 5 
With f=6? and using (17) and (19) we get 
51M U.)+M U,)]+0 6") 


©,” being defined in (18). 
This formula gives by simple means the value of 6,? correct to O (8°). 


(iv) The flow away from a stagnation point presents a formal difficulty in the 
integration of (16) which must be resolved, for if U=0O an infinite number of 
integral curves exists unless L(m)=0. Thus at a stagnation point we must have 
L(m)=0 or, from Table I, m= -—0.075. This value of m immediately fixes the 
value of @ at the stagnation point, and the value of 2 is found as follows. 

z 


Near m= - 0.075, L (m)=6 (m+0.075)+ a, (m+0.075% +... 


Hence from (16), when U =U, (x/c), 


2 
Lim (4%) 
dx 


z>0 U 
2 
=-6 Lim dx 
z->0 U’ 
~6 tim(#) 22 


+ The gradient of L(m) at m=—0.075 from Table I is not exactly 6, but is taken here as 6 
for convenience in fitting a later formula 
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TABLE If 
m H (m) L(m) 

+0.082 0 3.70 0.938 
+0.0818 0.011 3.69 0.953 
+ 0.0816 0.016 3.66 0.956 
+0.0812 0.024 3.63 0.962 
+0.0808 0.030 3.61 0.967 
+ 0.0804 0.035 3.59 0.969 
+0.080 0.039 3.58 0.971 
+0079 0.049 3.52 0.970 
+ 0.078 0.055 3.47 0.963 
+0.076 0.067 3.38 0.952 
+0.074 0.076 3.30 0.936 
+0.072 0.083 3.23 0.919 
+ 0.070 0.089 3.17 0.902 
+0.068 0.094 3.13 0.886 
. +0.064 0.104 3.05 0.854 
f + 0.060 0.113 2.99 0.825 
, +0.056 0.122 2.94 0.797 
+0.052 0.130 2.90 0.770 
+ 0.048 0.138 2.87 0.744 
+0.040 0.153 2.81 0.691 
+0.032 0.168 2.75 0.640 
+ 0.024 0.182 2.71 0.590 
+0.016 0.195 2.67 0.539 
+-0.008 0.208 2.64 0.490 
0 0.220 2.61 0.440 

— 0.016 0.244 2.55 0.342 
— 0.032 0.268 2.49 0.249 
— 0.048 0.291 2.44 0.156 
— 0.064 0.313 2.39 0.064 
— 0.080 0.333 2.34 — 0.028 
- 0.10 0.359 2.28 —0.138 
- 0.12 0.382 2.23 —0.251 
—0.14 0.404 2.18 — 0.362 
- 0.20 0.463 2.07 — 0.702 
- 0.25 0.500 2.00 — 1.000 


+ The solution of a laminar boundary layer is given on page 280 
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The values of 6? and its derivatives are therefore known and the integration 
of (16) can proceed from the stagnation point. 


(v) The distribution of /, m and H with x having been obtained, it may be 
required to construct the velocity distribution at any point. 


If the distribution is expressed in the usual way, u/U =f (y/8), it is very tedious 
to find expressions for the coefficients implied in this expression in terms 
of J, m and H. 


If however the distribution is expressed as 


the work becomes a great deal easier. This form of the velocity distribution is 
considered in detail in Ref. 14. It can be verified that: 


2 2 
 =m-- FO. and H= | F (dr, 
U y=0 y=0 


dy F’(0)’ \ay? [F’ (0) }° 
1 
l m 
whence F’ (0)= T° F” (0)= - H= F (t) dt. (20) 
Suppose that =F t<1. 
Then (20) gives at once a, = 2a, = 3 H 
2m _ 2 (4) 
whence 2 = + *) 


The velocity distribution can be calculated from this form, and the simplification 


resulting from the use of the form ; =F ( +) is clear. 


When a separation distribution is required, the form taken above for F (t) 
is insufficient, for then /=0 and the a’s take infinite values. This can be remedied 
by the use of the form F (t)=b,f'+5,t, for which it is easy to obtain 


(—2/m), b,=2H—(4/3)/(—2/m). 


An example of the use of this construction of distribution is given in a 
later paragraph. 
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Better approximations to the velocity distribution can be obtained by using 


3 
more of the boundary conditions (for example a) =0) and also by using 
y=0 


a function F which tends to infinity as (u/U)—> 1. However, since the approximate 
method gives only three pieces of information about the velocity distribution, it 
is doubtful whether better representations than those above could be obtained 
without very great complication. Even then, it is doubtful whether the approximate 
distribution would be sufficiently good for calculations such as stability calculations 
to be performed. 


(vi) It is worth while to remark here upon a test which is commonly applied 
to methods using the momentum equation but which has no significance, as will 
now be shown. In the flow U=U, (x/c}, a separation distribution is maintained 
for the value k= — 0.0904. 


The test applied to approximate methods is to assume a stream velocity 
U =U, (x/c¥ and find the value of k for which a separation distribution is maintained. 
The closeness of the value of k so obtained to - 0.0904 is taken as an indication 
of the accuracy of the method. 


However, the value of k so given depends only upon the value of H of the 
distribution which is assumed by the method to be a separation distribution. 


For suppose a distribution of constant “shape” is maintained for which 
~) =0, and so /=0. 
Oy] yao 
Then since the “shape” is constant, m has a constant value. Under these conditions, 
the momentum equation together with the boundary condition become, from (16): 


uu"! 


U d® d 
v (Uy 


On the assumption of U =U, (x/c¥, 


k (k -1)/k?=2 (H +2) or k= —1/(2H + 3). 


Thus the value of k obtained depends only upon the value of H assumed at 
separation, and clearly there is little significance in testing an approximate method 
in this way. 
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7. THE SIMPLEST METHOD 


(i) With the values assumed in Table I for /(m) and H (m), L(m) is found 
to be nearly a linear function of m, especially in the most common region of 7, 
say m>-0.1. This fact suggests a comparison of the L(m) curves for all the 
solutions known. Using the curves in Figs. 1-4 it is simple to calculate L (m) for ali 
the flows considered and in Fig. 5 the comparison is shown graphically. 


Two features of these curves are surprising. Firstly, it is curious that although 
the flows considered are of very different characters, the L (m) curves do not differ 
greatly from each other, and they differ by much less than the separate /(m) and 
H (m) curves, especially towards a point of separation. Secondly, it is curious that, 
at least in the range of m which is most common, i.e. m >-0.1, the curves are 
almost linear. 


These facts suggest strongly that a good approximation to all types of flow 
will be obtained by assuming a linear form for L (m) (or at most a quadratic form) 
in (16). By taking a reasonable mean of the curves in Fig. 5 we may assume 


It may be remarked that the possible range of choice for the two coefficients 
in (21), expressed as a percentage of their values, is about two per cent., which 
gives an indication of the error which may arise by the use of the form (21). 


A slightly more satisfactory form of L (m) would include the asymptotic suction 
distribution. For this, m= — 0.25, H=2.0, 1=0.5 and hence L (- 0.25)= -—1.0. The 
value of L (m) given by (21) at m= -0.25 is —1.05. A quadratic form of L (m) 
to satisfy the relations: L (-0.25)= —1.0, L (0)=0.455, L (0.075)=0.925, is 


L (m)=0.455 + 6.16m + 1.37m?. 


In the next section the consequences of the assumption (21) will be demon- 
strated, and similar consideration could be given to the quadratic assumption. The 
author feels that (21) is perfectly adequate. 


(ii) Substituting L (m)=0.45+6m in (16), we get 


which can be integrated to give 


0 
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—--—vu-=u, (2) 
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Thus (22) gives at once the distribution of 9 with x. At a stagnation point, it 
can easily be verified that (22) gives the values of 6? and (5-) which have already 
been given in Section 6 (iv). 

The distribution of m with x is known from that of 6° and then / and H can be 
calculated from Table I. Thus the local coefficient of skin-friction, 7, given by 


To 


pU* 


. i 
= and are obtained. 


Notice that although (22) was derived by a consideration of the momentum 


equation, that equation is not now required. The functional relationship between 
U dé U’?? 


leads at once to (22) and then only the / (m) and H (m) functions 


are required to obtain the skin-friction and displacement thickness. Now /(m) 
and H (m) assumed in Table I do not give a linear form for L (m) while (22) does, 
and so (22) in conjunction with Table I does not satisfy the momentum equation 
exactly. This is not necessarily a defect of the method, and indeed it may even 
give a greater “total” of accuracy on all r,, 6 and 6*. If for example a method 
gives inexact values for 7, and @, then satisfaction of the momentum equation will 
only lead to inexact values of 5*, whereas a different choice of H (m), while violating 
the momentum equation, might nevertheless give better values of 5*. We have 


therefore quite exhausted the possibilities of the momentum equation as a basis 
for approximate methods, and its prime and most important use must now remain 
as a test of the accuracy of exact and numerical solutions of the equation of motion. 

(iii) Young and Winterbottom) obtained a similar equation by making an 


approximation to, Pohlhausen’s method, their equation being (in our notation) 


x 


# =0.47vU aa | U*+***dx. It is clear from the foregoing analysis that this is 
0 
equivalent to choosing L (m)=0.47+6.28m. The gradient of this is rather too 


large, for when m > 0 this L (m) is rather greater than that for Schubauer’s ellipse, 
and when m <0 its value becomes rather too small. 


Falkner“) also derived an equation which in our notation becomes 
x 


2 
#=cv -* il U'dx | by varying the power of U occurring in the integral until a 
0 
best agreement was reached with the family of distributions given by U=U, (x/cy. 


It is not possible to derive an equivalent function L (m) from this, the accuracy of 
which it is difficult to estimate. 
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Tetervin’*) has given a general method of integration of the momentum 


\ay 


6 ‘ 
which Re= “ and k and n are constants. In the laminar flow case,-n=1 (see 


equation using the assumption that H is constant and 


in 


Section 7 (‘i)), and Tetervin’s assumptions are equivalent to assuming H (m) and 
1(m) to be constant, which is clearly insufficient for our purposes. From (15) it is 
clear that this assumption gives a linear form for L(m) so that the momentum 


equation is integrable to give a form for @ similar to (22). 


8. THE DIFFICULTY OF CONTINUOUS SUCTION 

This section is properly out of place in a paper which deals only with solid 
boundaries, but it is thought as well to explain the occurrence of a difficulty when 
there is a normal velocity at the boundary. Few exact solutions of flow exist with 
such continuous suction, and it certainly is desirable that an approximate method 
should be devised. An extremely brief discussion of the difficulties is given below. 


Putting y=0 in (1) and using the forms (8) we get 


Vv 


6v, 
— l(m), v, being the suction velocity. 


In Fig. 6 are drawn curves of +m-—(v,4/v)/1(m) against /(m), Table I being 
used in the calculation. (For the sake of simplicity in the argument, let us assume 
v, Varies inversely as 6 and that v.6/v= —0.3. The argument does not lose generality 
by this.) Each point on the curve @v,/v= —0.3 in Fig. 6 corresponds to a point on 
a boundary in a region of retarded flow in which separation takes place. Assume 
for simplicity that U’ is constant. Then for separation /(m)=0 and - 6*U’/v must 
equal +0.082. This is the first disadvantage of the present method, for separation 
is given by a constant value of —6°U’/v for all values of v., which cannot even be 
approximately true. Now as x increases, @ increases and also —U’6?/v. Thus x 
increases “up” the curve v.6/v= —0.3 in Fig. 6 from right to left. Eventually the 
maximum of this curve is reached, at / (m)=0.08, and — 6°U’/v must then decrease 
to reach the separation point. Thus the momentum equation must be so arranged 
that - 6?U’/v can decrease just at the right moment. Unfortunately this does not 
happen, and in any example a maximum value for —U’6?/v is reached beyond 
which the solution cannot be taken. Thus the separation value of U’6?/v is never 
reached and the method is unworkable. 


The real reason for the breakdown is the maximum in the curve of 
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m—(6v,/v)1(m) against /(m), which in turn is due to the zero gradient of the 
(m, |[(m)) curve at /(m)=0. For a method to be applied to continuous suction. 
the (m, 1(m)) curve must have a negative gradient at /(m)=0 and the solution will 
be valid only as long as 6v,/v is such that there is no maximum of m— (@v,/v) 1 (m) 
against /(m). Families of distributions which satisfy these conditions occur in the 


$(k-1) 
solution of U =U, (x/c)*, m <—0.0904 in which u varies as (*) and a paper 


is being prepared to show that these distributions do, in fact, resolve the difficulty. 


It is, however, very necessary to be able to estimate the accuracy of any 
approximate method, and for this reason an exact solution of a flow in which 
continuous suction is insufficient to prevent separation is urgently required. 
Watson"'*) has used his asymptotic theory to estimate the amount of suction 
required to produce a constant separation profile for the flow U=U, (x/c¥ but it 
is most necessary to have not only a solution for small velocities of suction, but 
also for some other flows. For example, a full solution of flow for U=8,-6,x 
for various values of v, would be of inestimable use. 


9. EXAMPLES 
Two examples will be given. 
(i) U=8,—,x, which has been solved by Howarth”). 


(22) gives immediately =0.075 ((1 - 1), é= 


Also, st 6? =m in this case, and /(m) and H (m) can be immediately evaluated 
from Tabk: I. Table II gives the various properties of the boundary layer for 


several positions of x. The whole work took less than an hour. 


IB, 1 /v 
6, 5* and \ B, 
against €. These quantities are non-dimensional. Also on Fig. 7 are shown the 
exact values obtained by Howarth. It is clear that the approximate method gives 
good results. The separation point given by m=0.082 is at €=0.117, whereas the 
exact result gives £=0.120. Thus the agreement is good, and’ could obviously be 
made better by assigning a slightly greater value to m at separation. This arbitrari- 
ness is at once a virtue and a vice of such an approximate method. 


In Fig. 7 are plotted the values of 


273 


| 


i 


| 
| 
| 
ajeuixoiddy yoexq 


a lo | 
| (uw) H ah | 


Il 


660°0 
161°0 
£07 


1780°0 
$980°0 
$8L0°0 
1990°0 
Lvro'0 
0L70'0 
£6100 
6$00°0 
0 


LITO 
O10 
$L80°0 
0$L0°0 
$790'0 
0s0'0 
0 


A 
oo 6 ¢ 6 6 = 
a @& 
nt 
Com NN TNH CO 
e a 
oeage 
222 
ack NA Ona No) 
| 
sta s 
o 6 6 6 6:85 6 6 
Nn 
Sis 
ol 
S 
wy 
274 


o-8 


07 


0-6 


os 


O-2 


0-24}— 


0-20}+— 


0-16 


0-08 


- 


APPROXIMATE CALCULATION OF THE LAMINAR BOUNDARY LAYER 


U= B.- Bi = | 


EXACT (HOWARTH) 


APPROXIMATE } 


° 
006 O08 O10 
Pe 


Fig. 7 


As a test of the accuracy of the construction of distributions given in (v) of 
Section 6, let us construct by that method the distribution at €=0.05 in this example. 
At this point, / (m)=0.178, m=0.0270 and H (m)=2.71. Hence from (20) a, =5.62, 
a, = — 2.39, a, =2.8 and the distribution can be quickly calculated. Table II tabulates 
some values, and also gives the exact distribution. Fig. 8 demonstrates the results 
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. . . 
and it is seen that the agreement is good, even though the approximate method gives f 
a finite value of y/@ at u/U=1. It is doubtful, however, whether this distribution ¢ 
is sufficiently accurate for the purpose of, for example, stability calculations. r 
(ii) Schubauer’s observed pressure distribution | 
The distribution of velocity over an ellipse which Schubauer observed has been i 
used as a test of many approximate methods. The validity of Schubauer’s results 
has been discussed by several authors and need nct concern us here since we are 
primarily interested in the accuracy of an approximate method of solving the 
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boundary layer equation. Hartree‘) has solved the equation on the differential 
analyser, using a distribution of velocity U which he obtained from Schubauer’s 


results; this may be regarded as an exact solution with which comparison is profitable. 


The complete solution according to the simplest method of this paper is shown 
2c 


2 
in Table Ill. The values of ) UU’ which are not denoted thus are 


2 
taken directly from Hartree’s tables. The values of ( é ) denoted thus + are 


5: UU’ and they in turn gave the values of U/U, 


denoted thus +. In the columns marked S, and S, are tabulated the summations 


obtained by integration of 


5 
occurring in the use of Simpson’s rule to integrate (7) . The approximate solution 


is started at x/c=0.2 since at that point Hartree started the analyser solution. A 
C2 
series solution was applied for < 0.2. The value of d *) is 
obtained from (22) and is given by 
0.2 


dx 


0 


5 
in which 6 and U take the values given by Hartree. Thus the integral | (7) a 
can be continued properly for x >0.2. From this U,@*/(vc) is obtained, and 
thence m, / (m) and H (m) from which can be directly 
obtained. The whole work as set out in Table III took almost one hour. For an 


arbitrarily given distribution of U /U,, whose gradient must be found, the work would 


take considerably longer. 


In Fig. 9, these approximate results are compared graphically with the exact 
results. The error in @ is less than 0.5 per cent. for x/c <1.75. As x/c increases 
from 1.75, the error increases to a final value of about two per cent. at x=1.95. This 
degree of accuracy in an approximate method is very good and justifies to a 
remarkable extent the use of the simple form (21). For x/c < 1.5, the error in 5* is 
less than 0.5 per cent., but for x/c > 1.5 the value of H (m) is systematically too 
high (as can be seen at once from Fig. 4) and thus in this range the error in 5* 
increases to a final value of about five per cent. The error in the skin friction, 


U, NU. oy y=0 
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Fig. 9 


is less than one per cent. for x/c < 1.4, but thereafter increases to a maximum of 
15 per cent. at x/c=1.8. At x/c=1.95 however the skin friction given by the 
approximate method takes the exact value. 


CONCLUSIONS 

A method has been derived in this paper for determining approximately, but 
with a good degree of accuracy, the principal characteristics of the laminar boundary 
layer. The method is extremely simple in use and has none of the disadvantages 
common to nearly all other approximate methods. 


Much of the paper discusses the uses of the momentum equation and a certain 
comparison of all known solutions of laminar boundary layer flow shows clearly 
the reasons for success or otherwise of an approximate method. It is also shown 
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that no advantage is gained by making more complicated initial assumptions when 
using the momentum equation, the possibilities of which as the basis of an 
approximate method are now exhausted. 


The solution of a laminar boundary layer (see Table I, page 265) is given by: 


x 


6? =0.45U U*dx . 
0 
U’?? 


Vv 


ou U 
= (m) 


= 0H (m) 
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